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I review the excursion set theory with particular attention toward apphcations to 
cold dark matter halo formation and growth, halo abundance, and halo clustering. 
After a brief introduction to notation and conventions, I begin by recounting the 
heuristic argument leading to the mass function of bound objects given by Press 
& Schechter. I then review the more formal derivation of the Press-Schechter halo 
mass function that makes use of excursion sets of the density field. The excursion 
set formalism is powerful and can be applied to numerous other problems. I review 
the excursion set formalism for describing both halo clustering and bias and the 
properties of void regions. As one of the most enduring legacies of the excursion 
set approach and one of its most common applications, I spend considerable time 
reviewing the excursion set theory of halo growth. This section of the review culmi- 
nates with the description of two Monte Carlo methods for generating ensembles of 
halo mass accretion histories. In the last section, I emphasize that the standard ex- 
cursion set approach is the result of several simplifying assumptions. Dropping these 
assumptions can lead to more faithful predictions and open excursion set theory to 
new applications. One such assumption is that the height of the barriers that define 
collapsed objects is a constant function of scale. I illustrate the implementation of 
the excursion set approach for barriers of arbitrary shape. One such application is 
the now well-known improvement of the excursion set mass function derived from 
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the "moving" barrier for ellipsoidal collapse. I also emphasize that the statement 
that halo accretion histories are independent of halo environment in the excursion 
set approach is not a general prediction of the theory. It is a simplifying assumption. 
I review the method for constructing correlated random walks of the density field in 
the more general case. I construct a simple toy model to illustrate that excursion set 
theory (with a constant barrier height) makes a simple and general prediction for the 
relation between halo accretion histories and the large-scale environments of halos: 
regions of high density preferentially contain late-forming halos and conversely for 
regions of low density. I conclude with a brief discussion of the importance of this 
prediction relative to recent numerical studies of the environmental dependence of 
halo properties. 



I. INTRODUCTION 



In the standard, hierarchical, cold dark matter (CDM) paradigm of cosmological struc- 
ture formation, galaxy formation begins with the gravitational collapse of overdense regions 
into bound, virialized halos of dark matter (DM). The average density of dark matter out- 
weighs that of baryonic matter by roughly six to one. Bound in the potential wells of dark 
matter halos, baryons proceed to cool, condense, and form galaxies. Understanding the 
fundamental properties and abundances of these dark matter halos is the first, necessary 
step in understanding the properties of galaxies. 

In this manuscript, I review the excursion set (also referred to as "extended Press- 
Schechter") approach to dark matter halo formation and halo clustering. I begin with 
some preliminaries that serve to define my notation and several common conventions in §1111 
I review the heuristic argument that Press & Schechter [l| used to derive their analytic 
mass function in § IIIII In this section I also draw attention to the most obvious weakness 
of the Press-Schechter argument, often referred to as the cloud-in-cloud problem. In § IIVI I 
review the theory of excursion sets of the density contrast field and the manner in which the 
excursion set theory solves the cloud-in-cloud problem. The standard excursion set theory 
then yields a halo mass function identical to the Press-Schechter proposal. 

The excursion set approach is very powerful and yields a great deal of insight into nu- 
merous aspects of halo formation and halo clustering. Following the mass function, one of 
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the most immediate applications of the excursion set approach is to predict the clustering 
properties of dark matter halos relative to the dark matter. I review the issue of halo bias in 
the context of the excursion set theory in § |Vl The bias relation derived in § |V] was derived 
prior to widespread use of excursion set theory, but this section also partially serves as a 
warm-up for the more complex arguments that follow. In § IVIl I summarize the excursion 
set theory of the void population. The content of § |V] and § IVIl contain some redundant 
information, but my intention is that repetition of the basic arguments in these sections will 
reinforce the logic and methods of the excursion set approach. 

In § IVIl I I review excursion set predictions for halo conditional mass functions, halo 
accretion rates, and halo formation times. The culmination of this review is the predominant 
algorithm for generating Monte Carlo merger histories on a halo-by-halo basis using the 
excursion set theory. Th,s algor.thm. refined most recently by Somerville & Kolatt (J and 
Cole et al. [3|] to similar effect, is the subject of § IVIIII 

Recent advancement in both theoretical methods and observational data have emphasized 
the weaknesses of the simplest implementation of the excursion set model. Nevertheless, the 
basic idea of the excursion set model is extremely valuable because it provides a toolbox for 
making quick estimates for numerous halo and galaxy properties. In addition, and perhaps 
more importantly, the excursion set approach remains one of the principle frameworks for 
qualitative reasoning and understanding of the complex results of direct, cosmo logical, nu- 
merical simulations. The algorithm for generating merger trees is one of the most widely 
used aspects of excursion set theory and is likely to remain one of the primary legacies of 
excursion set theory. These excursion set halo merger trees form the basis of countless semi- 
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orations of the physical processes involved in galaxy formation (examples are 



Refs. jsl, y, |5|). The excursion set halo merger trees are considerably less costly to produce 
than direct A^-body simulation, so they continue to be used to make general points (there 
are many examples, but an elegant recent example is Ref. j^) and to study the effects of 
varying cosmological parameters {e.g., Ref. Q])- Moreover, excursion set merger histories 
are also the basis of many models that aim to expand upon simulation results, either by 
extending results beyond the resolution of current simulations or building statistically-large 
halo samples in regimes where simulations have not, by modeling nonlinear dynamics in 

simplified ways {e.g., Refs. Q,y,y,Q,[iii,Q,Q)- 



The simplest implementation of the excursion set theory, which is the primary focus 
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of these lectures, is certainly not the only possible implementation of the excursion set 
approach. Excursion set theory can be modified to approximate numerous physical effects 
that are ignored in the derivations of the most basic set of equations. I address extensions of 
the excursion set models and proposed improvements in the final section of this review, § 11X1 
This section highlights areas where excursion set theory can be expanded and points toward 
the principle references where these extensions have been developed. I focus primarily on 
generalized forms for the barrier shape and excursion sets of correlated random walks. An 
interesting aspect of this section is the calculation of a general excursion set prediction for 
the dependence of halo formation times on the large-scale environments of halos. I close with 
a mention of the PINOCCHIO algorithm, developed by P. Monaco and collaborators primarily 

nnn 

in Refs. |14l . 1151 . 116I | , which makes use of Lagrangian perturbation theory together with some 
of the ideas of the excursion set approach to produce, of many things, halo merger histories. 

My primary aims are to present a relatively concise, pedagogical review of the basic logic 
and techniques of excursion set theory that can be used as a springboard to new studies and 
to collect the most important results of excursion set theory into a single manuscript, with 
a single notation. As such, the reference list contained in this manuscript is by no means 
an exhaustive review of the literature on analytic halo formation or even the excursion set 
approach. Many important contributions have been omitted in the interest of brevity and 
to limit the scope of the lectures. However, the references should be sufficient to develop the 
fundamental logic of the excursion set approach and the basic tools needed to extend this 
approach to new problems. As a minimal approach, the basic development of the excursion 
set approach can be followed through the series of papers by Bardeen et al. , Bond et 
al. 18 1, Lacey & Cole [lol,^, Sheth & Lemson |2l|, and Somerville & Kolatt j^. The various 
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manuscripts contained in the collection edited by Wax 
for understanding stochastic processes and deriving many of the basic results of excursion 
set theory (including results not contained in this review). Of the articles reprinted in the 



Wax collection, those by Chandrasekhar [23|, and Rice 2J, |25| are probably of the most 
immediate interest. 
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II. NOTATION AND CONVENTIONS 

In the following, I consider fluctuations in the density field p{x) described by the density 
contrast 5{x) = \p{x) — Pm]/pm, where pu is the mean mass density in the universe and x 
is a comoving spatial coordinate. In the standard paradigm, the universe is endowed with 
primordial density fluctuations during an epoch of cosmological inflation and the primordial 
density contrast is a statistically homogeneous and isotropic Gaussian random field. This 
means that the joint probability distribution of the density contrast at a set of points in 
space is given by a multivariate Gaussian distribution. Homogeneity requires that the mean 
{S{x)), of the distribution and the two-point function {5{xi)S{x2)) = ^{xi,X2) be invariant 
under translations. The two-point function is then a function only of the separation vector 
between two points, ^{xi,X2) = C,{xi — X2)- Isotropy requires that ^{x) is invariant under 
rotations as well, so the two-point correlation function is only a function of the distance 
between two points, ^{xi,X2) = Cd^i ~ ^2|)- 

The Fourier transform of the density contrast is given by the convention 

S{k) = J d^x (5(f)e^^-^ (1) 

with the inverse transform 

^^^^ ^ (2^ / "^'^ 5(^)e-^^^-^: (2) 

Notice that the 6{k) have dimensions of volume and that for a real- valued field 5{x), the 
Fourier coefficients obey the relation S{—k) = 6*{k). We have implicitly assumed that there 
is some very large cut-off scale L = V^l"^ that renders the integral / |5(x)|d^a; finite and that 
this scale is much larger than any other scale of interest so that it plays no meaningful role. 
Using these conventions, one can compute the two-point function ^(r) = -|- r)) in 

terms of the Fourier coefficients, where the average is taken over all space. The two-point 
function is a function only of the amplitude of r due to isotropy, and the result is 

i^r)-^,\k^V-^mf"^^^^^- (3) 
The correlation function is the Fourier transform of the power spectrum 



p{k)^v-wm?), 



(4) 



6 



where the average is over an ensemble of universes with the same statistical properties. The 
power spectrum has dimensions of volume and so a quantity that lends itself more easily to 
direct interpretation is the dimensionless combination 

A2(fc) = PP(fc)/27r2. (5) 

The correlation function is simply the mass variance. From Eq. ([3]), A^(/c) is the 

contribution to the mass variance from modes in a logarithmic interval in wavenumber, so 
that A'^{k) ~ 1 indicates order unity fluctuations in density on scales of order ~ k. 

In the standard, cold dark matter (CDM) model, A^(fc) increases with wavenumber (at 
least until some exceedingly small scale determined by the physics of the production of the 
CDM in the early universe), but we observe the density field smoothed with some resolution. 
Therefore, a quantity of physical interest is the density field smoothed on a particular scale 

6{x; Rw) = J d^x' W{\x' - x\; Rw)S{x') (6) 

The function W{x] -Rw) is the window function that weights the density field in a manner 
that is relevant for the particular application. According to the convention used in Eq. ([6]), 
the window function (sometimes called filter function) has units of inverse volume by di- 
mensional arguments. It is also useful to think of a window as having a particular window 
volume Vw- The window volume can be obtained operationally by normalizing W{x) such 
that it has a maximum value of unity and is dimensionless. Call this new dimensionless 
window function W'{x). The volume is given by integrating to give Vw = / d'^xW'{x). In 
this way, one thinks of the window weighting points in the space by different amounts. It 
should be clear that W{x) = W'{x)/Vw. Roughly speaking, the smoothed field is the av- 
erage of the density fluctuation in a region of volume Vw ~ -Rw- The Fourier transform of 
the smoothed field is 

6{k;Rw) = W{k;Ryv)S{k), (7) 

where W{k; -Rw) is the Fourier transform of the window function. 

The most natural choice of window function is probably a simple sphere in real space. 
The window function is then 

Wix;Rw)= ^ ^"-"""^^ (8) 
(x > -Rw) 
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In this case, the smoothed field is the average density in spheres of radius -Rw about point 
X. The window volume is simply Vw = 47ri?^/3. However, this choice of window has the 
undesirable property that the sharp transition in configuration space leads to power on all 
scales in Fourier space. Therefore, it is often convenient to smooth the boundary in real 
space. As I discuss below, it is often convenient to introduce particular window functions 
to ensure that the smoothed field has particular properties. The three most commonly used 
window functions are the real-space tophat window of Eq. ([8]), with Fourier transform 

Tj/^/ D N 3[sin(fci?i?w) - fc^wcos(A;/2w)] 

W{k; i?w) = , (9) 

the Fourier-space tophat window, and the Gaussian window. The Fourier-space tophat is 
defined in Fourier space as 

1 (k < Ry}) 

W{k-R^)= ^ - (10) 
{k > 

and is 

1 (sin(x/?w^) - xR;^^ cosjxR^^)) 

= 2^ ^''^ 

in real space. A disadvantage of this window is that it does not have a well-defined volume. 

This concern creeps up repeatedly in what follows. The Gaussian window is 

exp(-xVZR|0 /-.^N 

with a Fourier transform that also has the form of a Gaussian 

W{k- i?w) = exp(-A;V2i?w'), (13) 

with a width that is the reciprocal of its width in real space. The volume of the Gaussian 
window is Vw = {27t)^^'^R^. 

The density fiuctuation field is assumed to be a Gaussian random variable so the smoothed 
density fiuctuation field 6{x; R) is then a Gaussian random variable as well because it rep- 
resents a sum of Gaussian random variables. The variance of 6{x; R) can be computed in 
the same way as before [see Eq. ([3])] and is 

a\R) = {6\x;R)) = J dink A\k)\W{k; R)\^ (14) 

Thus, the probability of attaining a value of 5{x; R) between 6 and 5 + d(5 is 

P{6; R)d6 = , ^ exp[-6y2(r\R)]d6. (15) 
'2'Ka^{R) 
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FIG. 1: Power spectra in the standard ACDM cosmology with JIm = 1 — ^^A = 0.3, h = 0.7, 
(78 = 0.93, and rie/i^ = 0.022. The left panel shows the mass variance smoothed with a real space 
tophat window as a function of the smoothing mass or smoothing radius [Eq. ()14p ]. The right 
panel shows the rms density fluctuation per logarithmic interval of wavenumber as a function of 
wavenumber [Eq. ([3|)]. 



It is common to refer to the smoothing scale as either a length, as above, or a mass given 
by the mean density multiplied by the volume of the window M = PmVw. The left panel of 
Fig. [1] depicts a standard model power spectrum expressed as the variance of the smoothed 
density field, smoothed with a real-space tophat, and shows the smoothing scale in terms 
of both mass and length. The right panel of Fig. [1] shows the rms density fluctuation per 
logarithmic interval in wavenumber A{k). The power spectra were computed assuming 
= 1 — = 0.3, h = 0.7, cTg = 0.93, Q^h^ = 0.022, and a spectrum of primordial 
fluctuations with Pprim(^) oc A;. I have used the transfer function of Eisenstein & Hu 
to process the scale-invariant primordial spectrum into the post-recombination spectrum of 
density fluctuations out of which nonlinear structures form. 



III. THE PRESS-SCHECHTER MASS FUNCTION 

Press & Schechter [l| derived a relation for the mass spectrum of virialized objects from 
the hierarchical density fleld. In hierarchical models, there is structure on all scales and the 
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variance cr^(-R) — ^ oo as the smoothing scale i? — > 0. Press & Schechter essentially assumed 
that objects will collapse on some small scale once the smoothed density contrast on this 
scale exceeds some threshold value, but that the nonlinearities introduced by these virialized 
objects do not affect the collapse of overdense regions on much larger scales. The collapsed 
objects act only as resolution elements that trace the larger-scale fluctuations which may 
collapse at some later time. Strictly speaking, this is not correct; however, it is approximately 
true when primordial power spectra are sufficiently shallow that the additional large-scale 
power generated by nonlinearities is small compared to the primordial fluctuations on these 
scales (a study in one-dimension is Ref. 27|]). Moreover, this assumption leads to a simple 
parsing of the ingredients in the formation of nonlinear structure. The first ingredient is a 
characterization of the statistical properties of the primordial density fiuctuations. In the 
standard picture this is set during the infiationary epoch. The second ingredient is the 
evolution of overdensities according to linear perturbation theory. This is encapsulated in 
the growth function D{a) = 6(k, a)/S(k, a = 1) which is specified by the evolution of the 



background cosmology {e.g., Refs. |l9l . l28l . l29l|). The last ingredient is the threshold for 
collapse into a virialized object and is determined by examining the nonlinear collapse of 
spherical overdensities. 

Press & Schechter stated that the likelihood for collapse of objects of a specific size or 
mass {R oc M^/^) could be computed by examining the density fiuctuations on the desired 
scale. They continued by using a model for the collapse of a spherical tophat overdensity 
to argue that collapse on scale R should occur roughly when the smoothed density on that 
scale exceeds a critical value 6c, of order unity, independent of R. 

The implementation of the Press-Schechter prescription is simple. The mass within a 
region in which the smoothed density fiuctuation (dictated by the linear theory) is the 
critical value 6c, at some redshift z, corresponds to an object that has just virialized with 
mass M{R). The relationship between mass and smoothing scale is set by the volume of the 
window function. As an example, the relationship is M = 47rpM-R^/3 for a tophat window 
and M = {2n)^^'^pMR^ for a Gaussian window. Further, any region that exceeds the critical 
density fiuctuation threshold, will meet that threshold when smoothed on some larger scale 
R' > R. Consequently, the cumulative probability for a region to have a smoothed density 
above threshold gives the fractional volume occupied by virialized objects larger than the 



10 



smoothing scale F{M). Integrating Eq. (fT5|) . this probabihty is 



F(M) 



f mfl)d*=ierfc(^), (16) 

where erfc(a;) is the complementary error function, and u = 5c/<j{M) is the height of the 
threshold in units of the standard deviation of the smoothed density distribution. In this 
model, collapse of mass M is defined so that it occurs when the smoothed density fluctuation 
is 5c on the appropriate scale. Thus there is a typical scale that is collapsing at the present 
epoch, M^, when the variance is ct(M^) = 5c- 

In the hierarchical power spectra that we consider, (j{R) becomes arbitrarily large as R 
becomes arbitrarily small. Thus, -F(O) in Eq. (|T6l) should give the fraction of all mass in 
virialized objects; however, erfc(O) = 1 so that Eq. ( |T6l) states that only half of the mass 
density of the universe is contained in virialized objects. Press & Schechter noted this as 
a problem associated with not counting underdense regions in the integral Eq. (fT6!) . These 
authors argued that underdense regions will collapse onto overdense regions and multiplied 
F{M) in Eq. fll6p by a factor of two in order to account for all mass. Though the sense 
of this effect is certainly such that more mass will be contained in bound objects, that this 
should lead to precisely a factor of two increase in F{M) is far from convincing. 

Proceeding with this extra factor of two, the number of virialized objects with masses 
between M and M + dM is 



dM M 
In terms of the mass variance, this is 

dM = 



dF(M) 



dM 



dM. (17) 



dM 



2 pm Sc 


diner 


nM^ a 


dlnM 



exp(-|,) dM 

exp f - d^- (18) 
7rM2 dlnM ^ V 2 ; ^ ^ 

Without regard to the details of the shape of the power spectrum, cr(M) or z/(M), the mass 
function is close to a power law dn/dM oc M"^ for M <^ M^ and is exponentially cut-off 
for M> M,. 



IV. EXCURSION SET THEORY OF THE MASS FUNCTION 



A weakness of the Press-Schechter approach is that it does not account for the fact that 
at a particular smoothing scale 6{x; R) may be less than 6c, yet it may be larger than 6c 



11 



at some larger smoothing scale R' > R. It seems natural that this larger volume should 
collapse to form a virialized object, overwhelming the more diffuse patches within it. Clearly, 
the sense of including this effect is to increase the fraction of mass in collapsed objects and 
mitigate the factor of two discrepancy in the Press-Schechter formulas. In the literature, the 
issue of regions below threshold on a particular scale, but above threshold on a larger scale 
is referred to as the "cloud-in-cloud" problem. The cloud-in-cloud problem is closely tied to 
the Press-Schechter factor of two as I discuss below. 

To solve the cloud-in-cloud problem, it is necessary to compute the largest value of the 
smoothing scale for which the density threshold is exceeded. This was done in a formal way 
by Bond et al. [18|] (see also Refs. [30|, ImI, 1321]) and I follow their approach quite closely. The 
development of Bond et al. in turn follows closely the elegant review of stochastic processes 



by Chandrasekhar 



23| . and makes use of results from both Ref. 2J] and Ref. 33 1. 



In what follows, I consider the density contrast field evaluated at some early time far 
before any scales of interest have approached the nonlinear regime, but extrapolated to 
the present day using linear perturbation theory. In this way, the density contrast field 
does not need to obey the physical constraint 6 > —1 because this is merely the linear 
extrapolation of a density contrast of much smaller magnitude. In addition, all coordinates 
are Lagrangian coordinates defined at the same early time so that the position of each mass 
element is independent of time. The excursion set theory is, at its foundation, a set of 
rules for assigning mass elements to virialized objects of various sizes. Where needed, I 
will introduce the necessary conversions that map the initial Lagrangian sizes of overdense 
and underdense regions onto Eulerian coordinates, thereby accounting for the contraction 
or expansion of such regions respectively. 

Consider again evaluating the field 6{x; R) at various values of the smoothing scale R, at 
a single point x. For now, I will suppress the argument x and consider 6{R) as a function of 
smoothing scale at a single point in space. For very large R, (t{R) <^ Sc so the probability 
that the region lies above the boundary 6c is vanishingly small. With decreasing R, the 
standard deviation becomes larger and S{R) will eventually reach 6c at the first up-crossing 
of the boundary. The problem is to compute the probability that the first up-crossing of 
the barrier at 6c occurs on a scale R. For simplicity in what follows, let S = (J^{R) and 
let the value of S serve to denote the smoothing scale by exploiting the fact that S* is a 
monotonically decreasing function of R. I will then refer to the density contrast on that 
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scale as S{S). Throughout this discussion, it is important to remember that increasing 
S corresponds to decreasing R. The problem is to compute the probability that the first 
up-crossing through the barrier occurs between a value 5* and S + dS. 

Consider starting at a large smoothing scale, or small S = Si, where S{Si) = 6i < 6c- 
For a given change in the filtering scale AS*, there is some distribution for the probability of 
reaching 62 after an increment AS = S2 — Si > 0. In general, this probability distribution 
may depend not only on the size of the step AS, but the value of the density field on other 
scales. If the probability distribution for 62 after an increment AS" depends upon other points 
on the curve S{S), solving for the probability distribution of 5 at a given S is nontrivial. An 
important special case is when the smoothing window used to define S{S) is a k-space tophat 
as in Eqs. ( |TOl) -( fTTi) . In that case, increasing the filter scale corresponds to adding a set of 
independent Fourier modes to the smoothed density. These modes have not played a role in 
determining S{S) at other smoothing scales. In this special case, the transition probability 
for a change in density A6 associated with a change in filtering scale AS" is Gaussian with 
zero mean and variance S2 — Si = AS, independent of the starting point 5i(S'i). 



Following Bond et al. 



18(1, it is common to refer to a sequence of S{Si) given by many 



subsequent increases of the smoothing scale by increments A^i as a trajectory for S{S). In 
the case of k-space tophat filtering of the density field, each trajectory of S{S) executes a 
Brownian random walk. Three examples of such uncorrelated random walks are shown in 
Fig. [21 Notice that trajectories pierce the "barrier" at 6c many times and drop below 6c 
between subsequent up-crossings. The aim of the excursion set approach is to solve the 
cloud-in-cloud problem by determining the largest smoothing scale R or M, or equivalently 
the smallest value of the variance S, at which a trajectory penetrates the barrier at 6c. 

In the case of k-space tophat filtering, the probability of a transition from 61 to 62 = 
61 + A6 is 

n((52, ^2) d62 = ^iA6; AS) d(A5), (19) 

where 

is the Gaussian transition probability. Taking 5*1 = 0, 5*2 the smoothing scale of interest, 
and finding the probability of 62 > 6c returns the Press-Schechter probability for being in a 
collapsed object. The fact that some regions will exceed 6c for a smaller change in 5* and 



13 



CO 




S 

-M 



FIG. 2: Three examples of random walks of S{S) assuming each s step is independent as in the case 
of a sharp k-space window function. The axes are arbitrary. The horizontal dotted line represents 
some threshold value 6c- Notice that trajectories may penetrate the "barrier" at 6c many times. 



then fall below 6c by AS* has not yet been accounted for. 

Consider now, relating the distribution of 6 at one value of the smoothing scale 11(5, S) 
to the distribution on a subsequent step Il{6, S + AS) (smaller smoothing scale, larger S). 
This is 

U{S, S + AS) = J d{A6) ^(A(5; AS)U{6 - A6, S). (21) 

Taylor expanding Eq. (12T]) for small transitions, keeping terms up to (AS^), and integrating 
each term yields 

lim {A6)dU\ 
dS As^o \ 2 AS 05^ AS 05 )' ^ ' 

Using the fact that the transition probability is a Gaussian with {A5) = and ((A(5)^) = AS 

reveals 

^ = (23) 

dS 2 86^ ^ ^ 

as the relation governing the evolution of the probability distribution 11 with smoothing 
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scale. The probability distribution 11 for trajectories S{S) that never exceed Sc prior to 
some specific value of 5* can be obtained by solving Eq. fl23l) with the appropriate boundary 
conditions. The first boundary condition is that 11(5, S) is finite as 5 ^ —oo. Imagining 
each trajectory as absorbed and removed from the sample of trajectories at the moment it 
first crosses the threshold at 6c, then n(5c, =0 provides the second boundary condition. 

Eq. (!23l) may have a familiar form as it is also an equation describing diffusion with no 
drift or a one-dimensional heat equation describing heat flow in a long, thin, semi-infinite 
bar. The boundary condition at 6 = 6c is analogous to flxing the temperature to zero at 
this end of the bar. Take 6{So) = 6o as an arbitrary starting point for the random walk 
trajectories so that the initial condition is n(5o, 5*0) = 5d(5o) where 6b{x) is the Dirac delta 
function. Eq. ( 123|) can now be solved using familiar techniques. It is easiest to work in a 
shifted variable 7 = 5^ — 5 so that the flnite boundary is at 7 = 0. 

The first step is to Fourier transform Eq. (l23l) in the variable 7 (with conjugate u). Thus 
we define the Fourier transform of the probability distribution 

11(00, S) = Jd^ n(7, S)e'-\ (24) 

This gives an equation for the transform 



with solution 



U{uj, S) = c{uj) exp y-^Sj ■ (26) 
The boundary condition at 7 = {6 = 6c) guarantees that c{uj) is an odd function, so 

/•oo „ 

U{6, S) = / c{uj) sin(cu7) e"" (27) 
Jo 

Applying the initial condition to Eq. (p7|) gives 

c{uj) = ^ sin(u;7o) exp [y^o^ ' (28) 

where = 60 — 6c- The final solution then follows by inverting the transform, so 

2 ( — {S — So) \ 

11(7, S) = — sin(co'7o) sin(co'7) exp tu^ doo, (29) 

TV Jo \ 2 



TT JO 

from which 



/ (A5)n / [2i6c - 60) - A6f 
-^"P ^ 



(30) 
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where AS = S — Sq and AS = S — Sq. The solution to Eq. (l23i) in Eq. (130|) depends on the 
initial condition only with regard to the height of the barrier at 6c in relation to the starting 
position at 60 . The first term in Eq. (1301) is the Gaussian distribution that represents the 
points above threshold at S while the second term accounts for the trajectories that have 
been removed because they crossed above threshold at a filter scale S' < S but would have 
crossed back below the threshold by S. Chandrasekhar 23|| gives an elegant solution to this 
"absorbing barrier" problem almost upon inspection, exploiting the fact that random walks 
that hit the barrier have equal probability of continuing upward after their first encounter 
with the barrier as they do of continuing downward. Upon making that realization, the 
solution above can be obtained using the method of images and subtracting the contribution 
from another source beginning at 11(5, S = So) = 6b{S—26c), a shift of 6c above the threshold. 
The derivation given here is most useful because it contains the logic needed to generalize 
the problem to more complicated processes and barriers. 

The fraction of trajectories that have crossed above the threshold at or prior to some 
scale S{M) is the complement of the 11 distribution, 



F{S) = 1 - / n((5, S) dS = erfc 



( Sc-6o \ 



(31) 



\V2ASj' 

Taking 5*0 = and (5o = to indicate a starting value at very large smoothing scale, Eq. ( l3T|l 
is precisely the value given by Press & Schechter but without having to introduce an 
arbitrary factor of two to achieve the correct normalization. The second term in Eq. (1501) that 
accounts for trajectories that crossed threshold at some large scale and then traversed back 
back below threshold by S accounts for the probability missing from the Press-Schechter 
argument. 

The differential probability for a first piercing of the threshold then follows by differenti- 
ation, 



fiS\So,So)dS = ^dS 



- —00 

So 



dS 



{Sc - So 



Eq. (132|) is the fundamental relation as it gives the probability for first crossing a threshold 
given any starting point and any change in filtering scale A^*. The function / (5*150, 5*0) 
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is often referred to as the "first-crossing distribution." In cases where the random walk 
starts from 5*0 = with 6^ = these arguments are often suppressed and the first-crossing 
distribution is written simply f{S). The fraction of mass in collapsed objects in a narrow 
range of masses is obtained by setting So = and So = and rewriting Eq. ( 132|) in terms of 
the mass that corresponds to the variance M{S), 



dF 1 6r 



dM S 



dS 



dM 



2 



exp(-|l. (33) 



The mass function follows by using Eq. ffTTI) and substituting a = S results in 



dn 2 pm 6, 



dM V vr M2 a 



diner 



exp(-|j). (34) 



dlnM 

Again, this is precisely the Press-Schechter mass function with no ad hoc factor of two. 

Analytic mass spectra are often compared to simulation data without reference to specific 
power spectra or cosmological models by comparing the fraction of mass in collapsed objects 
per logarithmic interval in z/, vf{v) = dF/dlni/. For the standard excursion set theory, 

/M = /f-p(-y)- (35) 

In Fig. [31 I compare the Press-Schechter/Excursion set predictions for the fraction of mass 
in collapsed objects with the results of a suite of cosmological numerical simulations. While 
excursion set theory may explain the gross features of the mass spectrum of dark matter 
halos, it fails in its detail. The Press-Schechter/Excursion set relations predict too many 
low-mass halos and too few high-mass halos. Given the simplicity of the excursion set model 
described above, this is not surprising and the level of agreement provides encouragement 
that the excursion set model is a useful tool for understanding the gross features of halo 
abundance, formation, and clustering. In particular, these gross features are set by the 
statistics of the initial density fluctuations. 

One detail not discussed so far is the assignment of mass to a particular fllter radius R 
or variance 5*. As I have already mentioned, the Press-Schechter relation follows directly 
by assuming a window function that is a tophat in Fourier space, but such a window does 
not have a well-deflned volume to associate with it. The method of circumventing this that 
is implemented most often is to simply use the formulas derived from the sharp k-space 
filtering assumption, but assign mass according to a configuration space tophat filter so that 
M = 4:7iR^/3pM (though this is not the method used by Bond et al. in Ref. |18j). 
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0.01 




FIG. 3: Collapsed mass fractions. The solid line represents the standard excursion set theory 



M\ 



predictions. The dashed and dotted lines represent the improved fits of Sheth & Tormen 
and Jenkins et al. 35|] respectively. The points represent numerical data from a suite of N-body 



simulations reluctantly provided by J. L. Tinker. 



In the context of excursion set theory, the k-space tophat window function serves only 
to simplify the solution for 11(5, S) [Eq. fl30|) ] and the first barrier up-crossing distribution 
[Eq. (!32|) ]. In principle, other windows can be used but the mathematics become significantly 
more complicated. This is because the steps are no longer independent so it is necessary 
to compute the entire trajectory at once to account for the correlations between the steps. 
Ref. [l8] demonstrates the general procedure which amounts to generating a large number 
of solutions to a Langevin equation of random motion {e.g., , see Ref. j23|) with a corre- 
lated stochastic force. Each solution yields a single trajectory of density at each value of 
the smoothing scale. The probability distributions can be computed from an ensemble of 
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trajectories. This procedure is potentially complicated and time consuming, and it has the 
drawback that it does not yield a closed- form solution for the 11(5, S) or first up-crossing 
distributions. For these reasons, other filters are rarely discussed; however, it is important 
to keep in mind that the lack of correlations between different steps is not a prediction of 
excursion set theory. Contrarily, it is a simplifying assumption used in the most common 
implementation of excursion set theory. 



V. EXCURSION SET THEORY OF THE SPATIAL BIAS OF DARK MATTER 

HALOS 



It has long been appreciated that within the context of cosmological structure forma- 
tion theories, the clustering of dark matter halos differs from the overall clustering of mat- 



ter [34, 



36 



38 



391]. The excursion set formalism provides a neat framework with which 



to understand the relative clustering of halos. The argument is much the same as the 



peak-background split approach developed in Refs. |36l. l37l . |40|, and was expressed in the 



Jas I 



framework of the excursion set model by Mo & White [38|1 as I outline below. 

Consider the solution to the excursion set problem in Eq. (!30|) . This gives the probability 
distribution of 6 given that on a smoothing scale So, the smoothed density fluctuation is 
6o. Notice that the important quantity is the relative height of the density threshold 6c — So 
so that in regions with 5o > on large scales, trajectories are more likely to penetrate the 
barrier at 6c and conversely for 6o < 0. This means that Eq. fl3Tl) and Eq. fl32l) allow for the 
calculation of the mass function in regions with specific large-scale density fluctuations. 

The fraction of mass in collapsed halos of mass greater than M in a region that has a 
smoothed density fluctuation 6o on scale So (corresponding to mass Mq and volume Vq) is 
given by Eq. ( 13T|) . 

f6c 
; - 



F(M|<5o,^o)=erfc 



5n 



(36) 



V 2AS 

Notice that as the density of the region increases, F increases because smaller upward 
excursions are needed to cross the threshold. When 6o ^ Sc, F ^ 1 because the entire 
region will then be interpreted as a collapsed halo of mass Mq. The fraction of mass in halos 
with mass in the range M to M -|- dM is 

dS 



fiM\6o,So] 



dM 



dM 
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1 Sc - So 




exp 


[ {Sc - Sof 


A^3/2 


dM 


2AS 



dM, 



so that regions with smoothed density So on scale Sq contain, on average. 



J\fiM\So,So)dM=^ fiM\So,So] 



dS 



dM 



dM 



(37) 



(38) 



halos in this mass range. 

The quantity of interest is the relative over- abundance of halos in dense regions compared 
to the mean abundance of halos, 

AfiM\So,So) 



rL 

"halo 



{dn{M)/dM)Vo ^^^^ 
where dn(M)/dM is the mean number density of halos in a mass range of width dM about 
M from Eq. ( !34l) . The superscript L indicates that this is the overdensity in the initial 
Lagrangian space determined by the mass distribution at some very early time, ignoring the 
dynamical evolution of the overdense patch. 

The relative overdensity of halos in large overdense and underdense patches is easy to 
compute. In sufficiently large regions, So ^ S, Sq Sc- Expanding Eq. (!39|) to first order 
in the variables So/S and Sp/ Sc g ives a simple relation between halo abundance and dark 
matter density (see also Refs. jsJQ]) 



5, 



halo 



Sr 



-So, 



(40) 



where u = Sc/S^^'^ = Sc/(t{M) as before. The overdensity in the initial Lagrangian space is 

proportional to the dark matter overdensity and is a function of halo mass through u. The 

final ingredient needed to relate the abundance of halos to the matter density is a model 

for the dynamics that can map the initial Lagrangian volume to the final Eulerian space. 

Let V and S represent the Eulerian space variables corresponding to the Lagrangian space 

variables Vq and So- The final halo abundance is 

AfiM\So,So) , 
" {dniM)/dM)V ^ ^ 

Mo & White [ssl] give an extensive discussion of the mapping from Lagrangian to Eule- 
rian coordinates and use a spherical collapse model to determine the appropriate mapping 
quantitatively. In the hmit of a small overdensity So ^ 1, V Vo{l + S) , S So, and 



^halo = + 
= b^S. 



(42) 
(43) 



20 




FIG. 4: The large-scale bias of dark matter halos as a function of the scaled variable v. The 
solid line gives the prediction of the standard Press-Schechter or excursion set theory in Eq. ()42p . 
The dashed line is a modified form for the bias that follows from a modified form for the barrier 
criterion given by Sheth & Tormen in Ref. [34|. The dotted line represents an empirical fit to 



N-body simulation results given by Seljak & Warren 



39|. The points represent numerical data 



from a suite of N-body simulations kindly provided by J. L. Tinker. 



The halo overabundance is proportional to the matter overdensity and always has the same 
sense. This proportionality is commonly referred to has the halo bias 6h = (1 + [^'^ — l]/^c)- 
Interestingly, Eq. (1421) gives a mass scale at which halos are not biased with respect to 
the dark matter. At a halo mass equal to the typical collapse mass, M = M^, z/ = 1 or 
cr(M) = 6c so that the halo overdensity is the same as the total matter overdensity, 5haio = ^■ 
Larger halos will have larger values of u = 6c/cr{M) because a{M) decreases with mass, so 
larger halos cluster considerably more strongly than the overall clustering of mass, while 
halos with M < cluster more weakly than the overall mass distribution. 

The large-scale halo bias relation of Eq. fH2]) is plotted as a function of scaled peak height 
u, in Fig.m In that figure, I also plot numerical results for the relative clustering bias of halos 
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as well as improved fits for halo bias given by Sheth & Tormen 3^ and Seljak & Warren 39 1. 



The fit of Sheth & Tormen in Ref. 



34| can be motivated by excursion set theory with a 



modified form for the shape of the barrier as a function of smoothing scale {e.g., Refs. 41 



42|, 



and see § IIXI) . These simple prescriptions for halo bias form one of the cornerstones of the 
analytic halo model that is often used to quantify the clustering of both dark matter and 

E]. As 
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galaxies rapidly without resorting to a full numerical treatment [4c 
with the mass function, the standard excursion set approach gives an understanding of the 
general features of the halo bias, but fails to reproduce the halo bias in detail. 



VI. AN EXCURSION SET MODEL FOR VOIDS IN THE GALAXY 

DISTRIBUTION 



The excursion set model relates the initial density field on particular scales to specific 
structures in the evolved density field. To compute the properties of halos, one considers 
overdensities in the initial smoothed density field but it seems natural that this logic can 
be extended to underdense regions, and in particular, to voids in the density distribution. 
This is a particularly interesting application of the excursion set theory because of the 
relatively simple dynamics of underdense patches. Whereas overdense patches collapse, 
underdense patches expand more rapidly than the universal expansion rate. Contrary to 
the complicated behavior of overdensities, general underdense patches that are not initially 
spherically- symmet ric , topha t underdensities evolve toward spherically-symmetric tophat 
underdensities 
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5l| . This fact is encouraging. It suggests that a density threshold 



criterion based on a spherical, tophat model may be a more realistic model of void evolution 
than it is of the collapse of overdense patches. 

Sheth & Van Weygaert MJ made precisely this realization and applied the excursion set 
model to predict the abundance and sizes of voids in the universe. Most of the details of 
the excursion set model of voids can be found in their paper. The logic is quite similar 
to that of the previous two sections. First, it is necessary to assign a barrier height to 
designate a region that should be identified with a void, 6^ < 0. (At this point, I reiterate 
that the overdensities that we consider are initial overdensities evaluated at some very early 
time, before any scales of interest have gone nonlinear, and extrapolated to 2; = using 
the linear theory of evolution so that overdensities are not constrained to 6 > —1.) The 
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peculiar acceleration of a mass shell in a spherical underdense region is proportional to the 
integrated deficit in mass within the spherical region. A consequence of this is that inner 
shells accelerate outward faster than shells of mass on the outer edge of an underdense region. 
Mass shells start to accumulate in a narrow range of radii leading to a very high density 
ridge on the outer boundary of the underdense patch. Eventually, the inner shells overtake 
the outer shells at the moment of shell crossing. The relevant calculations are collected and 
reviewed in the appendices of Ref . jsil ■ 

In the now defunct standard cold dark matter model with Qm = 1, shell crossing for a 
spherical tophat overdensity occurs at a linearly-extrapolated overdensity of 6^ = —2.81. 
This is analogous to the case of spherical tophat overdensity collapse, where collapse to a 
point occurs when the linearly-extrapolated overdensity is 6c = 1.69. As with the spherical 
tophat overdensity, there is a weak dependence on cosmological parameters that I will ignore 
for the purposes of this discussion. Blumenthal et al. [5^ argued that galaxy voids should 
be identified with regions that have just experienced shell crossing. Again, this is analogous 
to defining halos as regions that have just experienced complete collapse. In this case, the 



galaxies in the void walls would 
model of Sheth & Van Weygaert 



ie on the overdense ridges that form at shell crossing. The 



5l| adopts this identification of voids with regions that just 



experience shell crossing and identify voids with the largest regions for which the smoothed 
density achieves the threshold 6^ = —2.81. 

In contrast to the case of halos, voids present an additional challenge to computing the 
first crossing distribution. The additional difficulty is that an underdense region may be 
contained in a still larger overdense region which may collapse around the underdensity, 
engulfing and eliminating the void. Therefore, it is sensible that there be a second criterion 
that the smoothed density not exceed some overdensity threshold 5^ > 0, on some scale 
larger than the first crossing of 6^. The appropriate height of the second barrier 6^ is not 
obvious. It is clear that 6^ < 6c because a void cannot reside within a virialized halo. Thus, 
the least restrictive choice that can be made is 6^ = 6c = 1.69. However, this does not 
account for the fact that some voids may be squeezed due to large-scale overdensities that 
have not yet collapsed and virialized and so it likely leads to an overestimate of the volume 
filling factor of voids in the universe. Sheth & Van Weygaert 5l| also explore 5^ = 1.06 
as this linearly-extrapolated overdensity corresponds to the moment of turnaround, or the 
transition from expansion to contraction, for a spherical overdensity. The drawback of this 
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choice is that it completely erases voids in regions that are only just turning around so that 
this choice likely leads to an underestimate of the filling factor of voids in the universe. This 
discussion suggests that the appropriate value of the "void-squashing" overdensity should 
be in the range 1.06 < < 1-69. 

The first step in modeling the void population with the excursion set theory should now 
be clear. It is necessary to solve for the first crossing distribution of the barrier 6^ subject to 
the additional restriction that the barrier at 5^ is not achieved at some larger scale either. 
As in § llVt it is convenient to solve for the probability distribution of trajectories that do 
not impact either barrier Il{6, S), and to derive the first-crossing distribution from Il{6, S). 
As before, the distribution of 6 at any value of 5* obeys 

^ = 1^. (44) 

dS 2 86^ ^ ^ 

In the case of voids, trajectories are removed from consideration when they impinge upon ei- 
ther of the two barriers at 5J or 6^ so that the boundary conditions are 11(6^, S) = n(5j, S) = 
0. 

This system of equations is the one-dimensional heat equation describing heat fiow in 
a long, thin bar of finite length held at fixed temperature at both ends. As such, it has 
a solution that is likely to be familiar. To begin with, unlike the case of the single-barrier 
problem in § IIVl the solution to this boundary value problem can only be a superposition of a 
discrete set of sinusoidal modes, namely those modes that individually obey the homogeneous 
boundary conditions. 

As before, it is convenient to define a new variable = S ~ S^. In that case, Eq. (jUj) still 
holds with 6 replaced by 7 and the boundary conditions are 11(7 = 0,5*) = 11(7 = 7c,v = 
5^ — dv, S) = 0. As before, trajectories start at 5 = at S* = so that the initial condition 
is 11(7, S' = 0) = (5d(7 + (5v)- The solution follows straightforwardly using separation of 
variables. Assuming that 11(7,5') = G{j)P{S), the functions G{-y) and P{S) obey 

1 dP _ 1 d^G _ 2 

PdS~2GW~ 

where k is some constant. Therefore P{S) oc exp(— /e^ S). The boundary conditions on G{'~f) 
can only be satisfied if > 0. Thus the solution is a superposition of sines and cosines. It 
is easy to show that the boundary conditions require that 

^(7) = !: c„ sin (—7). (46) 

^1 V7c,v J 
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The initial condition can be used to determine the c„ by using the orthogonahty relation for 
sinusoids, Jq sin(ny) sm{my)dy = {n/2)6nm, where 6nm is a Kronecker delta. The result is 
that Cn = 2sin(?7,7r|(5v|/7c,v)/7r. Consequently, the solution to the boundary value problem is 



5: - 5. 



E 

n=l L 



COS 



-6 — cos 



■[26,-5] 



X exp 



2(5- - 6, 



-S 



(47) 



In analogy with Eq. ( 13111 and Eq. ( l32l) . trajectories are removed from the region between 
the two boundaries at a "rate" of 

<5;f 



MS)dS 



IdU 



dS. 



(48) 



The lower limit of this divergence of trajectories indicates the trajectories that penetrate 
the lower boundary aX 6 = 6,. These are the void regions that we aim to model, so the first 
crossing distribution is 



MS)dS 



IdU 

y ^ —T^ — sm(?T,7rP) exp 

n=l 



n^TT^V^ S 



]dS, 



(49) 



where V = |5v|/(<5c + l^v|) is the void-in- cloud parameter, following the nomenclature of 
Ref. SlJ. The parameter V quantifies the relative importance of eliminating voids that lie 
within some larger overdense region in which the smoothed density contrast exceeds 6^. 

From the first-crossing distribution of Eq. (H9l) . it is straightforward to compute the 
"mass function" of voids, meaning the number of voids containing a specific amount of 
mass. As before, it is most convenient to represent this distribution in terms of a scaled 
variable = ISvl/S^^"^, so that the physical size distribution is determined by the size of 
the smoothing scale at a particular value of the variance 5* and the parameter T>. Then 
z/v/v(z/v) defined as the fraction of mass in voids per logarithmic interval dlnu, then follows 
by changing variables in Eq. fH9l) from S to as in § IIVI Explicitly, 



°° 2mxV'^ 
'^v/v('^v) = E — ~ — sin(n7rT') exp 

n=l 



2z/2 



(50) 



This void mass spectrum is shown in the left panel of Fig. [5] for the choices of 5, = —2.81 



and 6^ = 1.06 advocated by Sheth & Van Weygaert 



5l| . Two additional mass spectra 
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FIG. 5: Left: The mass distribution of voids predicted by the excursion set model of Sheth & Van 



Weygaert [5]|. The mass distribution is given in scaled units. The horizontal axis is the absolute 
value of the void threshold in units of the a{AI) for a particular mass scale z/y = \6v\/S^^'^. The 
vertical axis is the fraction of trajectories that first cross the threshold in a logarithmic interval of 
z^v Physically, this is the fraction of mass in voids that have a mass corresponding to a logarithmic 
interval in v^. The three lines show three choices of the overdensity threshold 5^. Right: The 
physical size distribution of voids. The three lines represent the number densities of voids of radius 
R for the same three choices of 6^ as in the left panel. 



derived from different values of 6^ are shown in Fig.[5]to illustrate the importance of the void- 
in-cloud parameter. The choice 5^ — oo corresponds to ignoring completely the possibility 
for a larger-scale overdensity to overwhelm a smaller void within it. Decreasing 6^ makes 
the probability for a larger overdensity to overwhelm a void within it higher and decreases 
the number of voids of all sizes. Two things are worthy of note. First, the number of large 
voids (high u^) is independent of This is because it is very unlikely for a trajectory to 
cross both the 5^ and the 6^ thresholds at small values of S. Second, regardless of the choice 
of 6^, the fraction of mass in voids is peaked around voids with z/y ~ 1 — 1.5. This places 
the typical void mass in the broad range M^oid ~ 10^^ — 10^^ /i^^Mq. Of course, as with 
the mass function, the number of voids of a particular size will be shifted toward smaller 
objects because a fixed fraction of mass can be partitioned into more small voids and fewer 
large voids. 
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The void size spectrum in Eq. ( !50|) gives the fraction of mass, or of volume in the initial 
Lagrangian space, in voids of a particular size. The size variable is the scaled variable u^. 
To convert this to a physical size spectrum where the sizes of objects are measured in units 
of length or volume requires relating the initial Lagrangian volume to an Eulerian volume. 
The simplest conversion follows from the fact that in the spherical tophat model, the tophat 
radius expands by a factor or ~ 1.72 by the time shell crossing occurs. Finally, one also 
needs to assume a power spectrum to relate particular values of z/^ with volume in the initial 
Lagrangian coordinate system. As before, I will assume the spectrum in Fig. [H The logic 
is then the same as that of the mass function. The number density of voids as a function 
of physical void radius is shown in the right panel of Fig. [5] and reflects the basic features 
of the voids model. The number density of voids is independent of 5^ for large voids, and 
the number of small voids is suppressed because small, underdense regions are likely to be 
contained within larger overdense regions. 

So far, I have ignored the primary shortfall of the void model. The weakness of the 
approach of Ref. 511] is that it defines a void as a region that is underdense with respect 
to total mass whereas observationally voids are identified as under densities in the galaxy 



distribution, as noted in Ref. 



5l|. Furlanetto & Piran 



53| have extended the Sheth & Van 



Weygaert model to identify voids as regions that are underdense in the galaxy distribution. 
To do this, these authors added to the excursion set theory a method for assigning galaxies 
to halos and thereby established a mapping between the galaxy number density in a region 
and a linearly-extrapolated overdensity. Methods for assignin g ga laxies to halos are beyond 
the scope of this article so I will not discuss the model of Ref. 53|] any further except to say 
that it is a more practical approach to voids because it phrases predictions in terms of the 
observed galaxy overdensity. 



VII. HALO FORMATION IN THE EXCURSION SET THEORY 

Excursion set theory is also a powerful formalism with which to study the formation 
history of halos. It is most natural to picture evolution occuring through the linear growth 
of the density field. In that way, the smoothed overdensity scales with time according to 
the linear growth factor D{a). Writing the smoothed density field as a function of both 
the smoothing scale S and the cosmological expansion factor a, this evolution is 6{S, a) = 
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D{a)S{S,a = 1). As a matter of convention, D{a) is normalized so that D{a = 1) = 1. 
In the framework of excursion sets, predictions depend upon the density field in the ratio 
u = 6c/(y{M) and it is often easier to envision the density field as fixed at S{S) = 6{S, a = 1) 
and to consider the height of the critical density threshold as a function of time. In that way, 
collapse at a scale factor a' 7^ 1 corresponds to the a = 1 density fluctuation penetrating a 
barrier of height Sc{a') = 5^/ D{a'). 

In their pioneering paper. Bond et al. jisj] briefly discussed the utility of the excursion 
set a ppr oach for understanding the formation histories of halos. Subsequently, Lacey & 



Cole [19|] studied this problem in great detail in a rather elegant paper and I larg ely follow 
their treatment in what follows. I will adopt the simplifying notation of Ref. 19|] and refer 
to the barrier height for collapse at a specific time as 

uj{a) = 5jD{a). (51) 

The variable uj indicates the time of the collapse through only the linear growth function 
D{a). In the "old standard" CDM model with VLm = 1, the growth function is simply 
D{a) = a and the conversion between 00 and a is trivial. However, in the currently-favored 
ACDM model, the growth function is given by an incomplete beta function integral 2^ and 
the conversion is more complicated. 

Many of the insights into halo formation can be gleaned by considering the two-barrier 
problem, just as with halo bias and void regions in the previous two sections. Again following 



the work of Lacey & Cole Jjj], I also adopt a simplifying notation for the conditional 
probabilities of barrier crossings. I illustrate this notation by re-writing the two-barrier 
probability distribution already given in Eq. fl32l) and Eq. fl37|) as 



f{S2,UJ2\Sl,UJl)dS2 = ^^^/2 6XP 



2AS 



dS2, (52) 



where Au = 002 — 001 is the difference between the two barrier heights (assumed to be 
arbitrary with 002 > ooi), AS = S2 — Si. As always, some mass Mi and M2 are associated 
with the variances Si and S2 respectively. The quantity /(S'2, oo2\Si, co'i)dS'2 is the conditional 
probability that a trajectory pierces the barrier 002 in an interval of width dS'2 about S'2 on 
the condition that the trajectory first pierces ooi at 5*1. As a further shorthand, I take 
/('S'2, ^^2) dS'2 to denote the probability starting from a very large smoothing scale Si = 0, 
and ooi = 0. 
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FIG. 6: The conditional mass function. All panels refer to the conditional mass functions of halos 
fixed at z = looking back toward higher redshifts as labeled in each panel. Four halo masses of 
interest are shown. 

In this section and those that follow, the interpretation of Eq. (1521) is different from that 
in § IIVI and § |Vl In those sections, the initial condition for the random walk was fixed by 
the large scale environment on some larger filtering scale Si ^ 5*2. Here, the two barriers 
represent the critical density for collapse at two different times so the change in barrier 
height Allj represents a shift in time. 

A. The Conditional Mass Function 

The most immediate quantity of interest is the conditional mass function. Given a halo of 
mass Ml at time ti, one can compute the average manner in which this mass was partitioned 
among smaller halos at some earlier time t2 < ti {002 > "-^i)- The conditional mass function 
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is simply the average number of halos of mass M2 at time ^2 that are incorporated into a 
mass Ml object at time ti. This is precisely the two-barrier problem. A halo of mass Mi 
has its mass partitioned, on average, among a spectrum of halos at time W2 given almost 
directly from Eq. ( l52l) {e.g., see Refs. 

■ d& 



dn(M2|Mi) Ml 

- /(S2,cc;2|Si,tx;i) 



dMo 



dM2. (53) 



dM2 M2 

The function /(S'2, cc'2|S'i, cui) gives the probability of the second barrier crossing at a par- 
ticular value of 5*2, while the factor M1/M2 converts this from a probability per unit mass 
of halo Ml into the number of halos of mass M2 . 

Figure [6] shows conditional mass functions for standard CDM power spectra normalized 
such that M^ = 10^'^ /i~^Mq (erg — 0.93) for four halo masses at four different lookback 
redshifts, z = 1/4, z = 1/2, z = 1, and z = 3 (all panels have final redshift at 2; = 0). The 
evolution toward increased fragmentation as the second epoch moves to higher redshifts 
is evident. Notice that more massive halos fragment into small sub-units more quickly. 
Also, notice that each conditional mass function goes toward an approximate power law as 
M2/M1 -C 1. This can be understood directly from Eq. ([6]). Assume that the power spectrum 
is an approximate power law S{M) ~ 6'^{M/M^,)~^ in the mass range of interest. In the limit 
of a large mass difference, Auj/AS < 1 and M2/M1 < 1 and dn/dM2 oc M2^^~^^^^^. Near 
M^, CDM-like power spectra vary slowly, with an effective power-law index near /3 ^ 0.35 
giving dn/dM2 oc Mg^^^^, which is approximately the power law shown in Figure [61 

B. Halo Accretion Rates 

The excursion set theory also provides a framework for computing mass accretion rates 
as well. The two-barrier problem of § IVII Al can be manipulated to yield an average mass 
accretion rate. Consider the probability of transitioning forward in time, from M2 at barrier 
002 to Ml at uoi. Bayes' rule for combining probabilities gives the forward probability as 13] 

f/C IC \AC f{S2,UJ2\Sl,UJl)f{Sl,UJi) 

f{Si,uJi\S2,uj2)dSi = — r d5^i (54) 

J[02, ^2) 



1 cui(u;2 - u;i) ^ "'^'^^ 



27r UJ2 



S2 

Si{S2 — Si 



(55) 



/ {UJ1S2 ~UJ2Si) \ 



30 




FIG. 7: Halo accretion rates per logarithmic interval of mass change. This is based on a very 
similar plot in Ref. 



19l |. The plots show the accretion rates per unit mass change for halos of 
five different masses (shown in the left panel) at z = 0. The left panel shows the probability for 
a change AM as a function of AM [Eq. ()58p ]. while the right panel shows the fractional mass 



accretion rate given by (AM/M2) x d^i?/dln AMdlna. 



In language familiar to cosmology, f{Si,uJi\S2,uj2) is the posterior probability, f{Si,uJi) is 
the prior, and f{S2,uj2\Si,uJi) is the likelihood. 

The instantaneous probability per unit shift in barrier height follows by taking Au = 
{uj2 — cui) — > in Eq. ( l5il) . giving 

-1 3/2 



duj dSi 



So 



exp 



2S2S1 



(57) 



SiiS2-S,) 

from which the probability per unit change in mass AM = Mi — M2, per unit time (or 
dlna = H{a)dt) is 



dlnAMd In a 



'2 AM Lj/a{Mi) 

TT Ml (1 - Si/S2)3/2 



X exp 



coHS2-Si) 



d Incj 



dlna 



diner 



din Ml 



(58) 



2^1^2 

The average rate of mass accretion per unit change in mass can be obtained by multiply- 
ing Eq. ( l58i) by AM, while the average total mass accretion rate follows from multiplying 
Eq. (!58l) by AM and integrating over In AM. 
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FIG. 8: Total halo mass accretion rates in units of halo mass and as a function of halo mass at 
four different redshifts. 

The excursion set accretion rates of halos in the standard ACDM cosmology are sum- 
marized in Fig. [7] and Fig. [HI The left panel Fig. [7| shows the probability for a particular 
fractional mass change as a function of the fractional mass change. The rate of accretion 
by small mass increments diverges. However, the average rate of mass increase (in the right 
panel of Fig. [7]) converges and is dominated by the less frequent, high- mass mergers. Notice 
in particular that low-mass halos increase their masses very quickly via mergers with halos 
larger then themselves while halos with experience only infrequent mergers with 

halos larger than themselves. 

The total mass accretion rate given by integrating Eq. ( l58l) over In AM is shown in Fig. [8] 
for various redshifts as a function of halo masses. Fig. [8] shows explicitly that low-mass 
halos should be expected increase there mass many times over in a single Hubble time. The 
standard interpretation of the collection of results shown in Fig. [7] and Fig. [8] is that low- 
mass halos (M < M^) should be expected to increase their masses by mergers with halos 
larger than themselves, while high-mass halos (M > M^) typically absorb smaller halos and 
increase their masses relatively slowly. 
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C. Halo Formation Times 

It is often useful to have a sense of a time by which a halo acquired most of its mass. For 
example, one may be modeling galaxies and may try to estimate stellar evolution histories 
within a particular halo. Several definitions of halo formation time exist. An intuitive 
definition is the time that a halo first acquires half of its final mass because once a progenitor 
of a halo has achieved this threshold it can uniquely be identified as the main progenitor 
of the final object. Consider the formation time given by this definition uj^ or, in terms of 
expansion factor, ap. 

It is not trivial to predict halo formation times from excursion set theory by considering 
the two-barrier problem as before. The reason for this difficulty is simple to understand. 
Consider the formation of a halo of mass Mi at time Ui. The probability of crossing a 
higher threshold UJ2 (corresponding to earlier time) at some value of the variance greater 
than 5*2 = S{M2 = Mi/2) > Si gives only the probability per unit mass that the halo 
has some progenitor with mass M2 < Mi/2 at time 002- It does not guarantee that this is 
the mass of the most massive of all of the subunits that merged to form the halo of mass 
Ml at uJi. Lacey & Cole 191] circumvent this difficulty by giving an approximate counting 



argument for the formation time. I reproduce the argument of Ref. [19| in what follows. 

First, the results of § I VI I Bl give the number of halos of mass M2 at time ^2 that are 
incorporated into halos of mass Mi at ti > t2 as 

d^n dn(M2) 
dM^r^MT^^^'^'"'^^'^^'^- ^''^ 
So long as M2 > Mi/2, each trajectory must connect two unique halos because there cannot 

be two paths each of which contain more than half of the final halo mass. However, it is 
possible that a halo of mass Mi at ti has no progenitor of mass greater than Mi/2 at time 
^2- The counting argument accounts for this. The probability that a halo of mass Mi at ti 
has a progenitor in the mass range Mi > M2 > Mi/2 at time t2 is then given by the ratio 
of halos that evolve to a halo of mass Mi relative to the total number of halos of mass Mi, 

dP 



[dn{M2)/dM2]fiSi,ui 


S2, 




[dn{Mi)/dMi] 


|dSi/dMi 


-1 



Using Bayes' rule again and the definitions of the mass functions in terms of the first-crossing 
distributions gives 



dP _ /Ml 

dMs ~ VM2 



f{S2,UJ2\Si,UJi] 



dS2 
dMo 



(61) 
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FIG. 9: The probability distributions for formation times in the standard ACDM cosmology with 
a scale-invariant primordial power spectrum normalized to ag = 0.93, so that = 10^^ /i^^M©. 
The left panel shows the distribution as a function of threshold height dP{uJF)/dujp while the right 
panel shows the probability distribution of formation scale factors dP(oF)/daF. 



The cumulative probability of a formation time prior to time t2 is then the integral 

rS2=S{AIi/2) /Mi\ 

P(tF < ^2) = P{0OF > UJ2) = J^^ ij^)f{S',,UJ2\SuUjMS2, (62) 

and the differential probability can be obtained by differentiating Eq. (!62|) with respect to 
either up or tp. 

Formation time distributions computed from Eq. (1621) are shown in Fig. M for a variety 
of halo masses at 2; = (a = 1, tup = (^c — 1.69). The most obvious feature of formation 
times is that small halos form early and large halos form relatively late in accordance with 
our general expectation in hierarchical models of structure formation. Note also that in the 
simplest implementation of the excursion set formalism, in particular making use of k-space 
tophat filtering of the density field, the conditional probability /(S'2, co'2|S'i, cJi) is a function 
only of the differences S2 — Si and 002 — uji- In this case, this means that the formation times 
in particular and, more generally, the entire formation histories of halos are independent 
of the larger-scale environments that halos reside in. In the next section, I discuss halo 
formation times in the context of halo merger trees. 
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VIII. HALO MERGER TREES 

The probability of first up-crossing at a second barrier given a particular starting point is 
given by the solution to the two-barrier problem in Eq. ( l52l) . This relation has been repeated 
several times throughout these notes in several forms because it is the fundamental relation 
on which most physical results are based. In the preceding sections, I have applied this 
formula to derive relations about various aspects of halo fragmentation and formation over 



time. So far, these quantities have been averages 



as in § IVII Al and § IVIIB|) or probability 



distributions (as in § IVII CI) . However, as Ref. 19|] describes, the interpretation of Eq. ( !52|l 
as a probability for a transition from one barrier to another means that one should be able 
to draw specific fragmentation probabilities from this distribution repeatedly for a number 
of time intervals Acjj, and thereby reconstruct individual trajectories of halo mass versus 
cosmic expansion factor. This enables one to follow the fragmentation of halos on an object- 
by-object basis. This is the logic behind the algorithm for the generation of Monte Carlo 



merger trees described by (see also earlier methods in Refs. 54, l55|, l56l|). 

The trajectories 5(5") have structure on arbitrarily small scales corresponding to mergers 
with very small halos. This is a restatement of the divergence of the mean number of 
transitions with small mass changes shown in the left panel of Fig. [71 Consequently, it is 
necessary to examine trajectories with a particular resolution in order to "smooth over" 
the mergers with numerous tiny halos. In practice, each application to the prediction of 
a specific physical quantity has a minimum mass scale of interest Mmin, so it is natural to 
set the resolution with which trajectories are computed according to this minimum scale. 



Consider Eq. ((521) in the limit Acu/VAS < 1. This gives a transition probability / 
Aco'/(\/27rAS'^/^), that is directly proportional to the time interval Acij. Lacey & Cole [l 
interpreted this as indicative of a probability that is due to a single, improbable merger event. 



For this reason, Lacey & Cole refer to the limit Auj/^/ AS ^ 1 as the "binary regime." This 
can be used to motivate a particular choice of step-size Auj that sets the resolution of the 
trajectory. 



First, consider making the change of variable to m = Acj/V AS in Eq. (1521) . This converts 
the probability distribution to a Gaussian in the variable u with a mean of zero and a unit 
variance. Of course, the probability distribution is one-sided in that it is subject to the 
constraint u > 0. To resolve individual, binary mergers of a parent halo of mass M with a 
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merging halo at the threshold M^^^ the step-size should obey 

M^i„ (63) 

so that these mergers are in the binary regime. Authors often parameterize this guideline 
and set Atu = /step \J\dS/dM\ Mmin • The guideline in Eq. ( 163!) makes a transition with 
AM > Mmin a rare event. The motivation behind this choice is that it is very unlikely that 
the transition could be due to two rare events, so that most mergers with objects of mass 
> Mmin are resolved as binary encounters. In particular, the probability for a transition 
with AM > Mmin is 

2 /"/stop/v^ 

P(> AM) = / exp( V)dy 

VTT JO 

= erf(/,tep/v^) (64) 

/step ) 

where the last step assumes /step <^ 1- The choice of /step is motivated by the need to resolve 
all mergers of objects of mass Mmin, driving /step toward small values, and the desire to spend 
as little effort as possible computing the probabilities for transitions with AM < Mmin, 
driving /step toward large values. As I discuss below Ref. advocates /step~ 0.3 in order 
to gu^antee accuracy down to M„,„. TVpica, value. a.e /,..„<ai. 

With these preliminaries out of the way, the Lacey & Cole [19| algorithm for generating 
merger histories is as follows. First, determine the appropriate timestep using Eq. fl63l) . 
Second, select a transition AS* from the probability distribution of Eq. (!52|) and invert the 
S{M) relation to obtain both the change in mass AM and the new main progenitor mass 
M' = M — AM. One then repeats this procedure at the new values of u and S{M') in order 
to determine the next fragmentation of the main progenitor. This process continues until 
the remaining mass is less than Mmin- 

The trajectory for the main progenitor, defined as the most massive progenitor at each 
timestep according to this algorithm, may form the trunk of a "merger tree." All of the 
merging halos of mass AM have prior fragmentation histories of their own. For each mass 
AM above the threshold Mmin, one can generate an independent history for the infalling 
"branch" on the tree using the same algorithm. 

The process of constructing a mass accretion history or merger tree in this way must be 
repeated numerous times in order to sample the variety of ways in which a halo of a fixed 
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mass at a fixed time might build up its mass. This is a Monte Carlo method for exploring the 
various mass accretion histories. It is conventional to refer to each individual tree generated 
in this way as a particular realization in an ensemble of merger histories. 

The algorithm for generating Monte Carlo merger trees described above is convenient 
because of its simplicity. Unfortunately, while this algorithm conserves mass by construc- 
tion, it overpredicts the number of progenitor halos with masses M > M^[^ at previous 
timesteps relative to the analytic distribution of Eq. fl53l) . This point has been emphasized 



by Somerville & Kolatt (see also Ref. 21|). In other words, the Monte Carlo procedure 
does not lead to a mean population of halos at high redshift that is consistent with the 
excursion set relations of the previous section. In the left panel of Figure [lOl I reproduce 
Figure 2 of Ref. 21] in order to emphasize this point. The smooth lines in the left panel 
of this figure represent the analytic conditional mass functions described in § IVII Al and 
the histograms represent conditional mass functions from an ensemble of merger histories 
constructed according to the method of Lacey & Cole 13]. That the merger tree method 
of Lacey & Cole becomes more and more inconsistent with the analytic conditional mass 
function with increasing redshift is evident. 

Somerville & Kolatt [2] suggested two very plausible sources for this discrepancy. The 
first is that after choosing the first mass AM from Eq. (!52|) . the mass of the second body in 
the merger is forced to be M' = M — AM. This choice is made in order to conserve mass, 
but without regard to the probability of a transition to the mass M' . This drives an over- 
abundance of high-mass halos because all of the remaining mass is forced into a halo of mass 
M'. The second source stems from the fact that the simple algorithm above largely neglects 
accretion with transitions AM < Mmin when the probability for such a transition becomes 
infinitely large as AM — > 0. While Ref. ^ does not present a rigorous demonstration of a 
self-consistent algorithm, they give a qualitatively well-motivated algorithm that reproduces 
analytic halo conditional mass functions at high redshift very well. This new algorithm is 



now the most commonly used method and 
is, in practice, very similar to that of Ref. 



Sheth & Pitman 



summarize it below. A successful algorithm that 
is developed through the papers of Sheth js^] , 



58| . and Sheth & Lemson 



2l| . though as far as I am aware, this approach 



also lacks rigorous support other than that it is exact for clumps in an uncorrelated density 
field. Cole et al. [3| provide another variation of the Lacey & Cole |l9] merger tree algorithm 
that is similar in spirit and effect to the that of Ref. 



37 






Z = 1.0 

i 1 1 'i^'- 


i 1 1 "III 

;\. z = 3-0 I 

1 li ... 


z = 6.0 ; 

i\ ^ 



10 100 



FIG. 10: These figures are reproduced from Ref. [4]. The left panels show the distribution of 
progenitor masses logdA^/dM for a halo of initial mass Mj = 500 at four different redshifts. The 

n 

histograms were computed using the binary merger tree method of Lacey & Cole in Ref. [19|] and 
the smooth lines are the conditional mass functions of excursion set theory. The mass units are 
arbitrary. The discrepancy between the conditional mass functions and the Monte Carlo merger 
tree results is clear and becomes larger with increasing redshift. The right panels show the same 
physical quantities. In this case the smooth lines are the same conditional mass functions as before. 
The histograms show the results of the improved Monte Carlo merger tree algorithm proposed by 
Somerville & Kolatt [2]. The discrepancy between the Monte Carlo method and the analytic 
conditional mass functions is clearly reduced using the Somerville & Kolatt modification to the 
merger tree algorithm. 



The Somerville & Kolatt 121 alErorithm beerins in the same manner as the original Lacey & 
Cole [191 procedure. First, choose a timestep corresponding to the desired mass resolution in 
accord with Eq. fl63l) . At each timestep choose a transition AS* from Eq. (l52l) and convert it 
to a mass change AM. The algorithm now deviates from that of Lacey & Cole because the 
Somerville & Kolatt algorithm allows for both multiple mergers and diffuse mass accretion. 
At this point, I note that the term "diffuse mass" is used simply to refer to accretion of 
objects with mass below the threshold Mmin- In the standard excursion set approach, there 
is no truly diffuse mass as all mass is bound into halos of some size. 

The Somerville & Kolatt algorithm continues as follows. If AM < M^in, remove it 
from further consideration and regard it as diffuse accretion. If AM > Mmin treat it as a 
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progenitor halo that has merged. Next, update the remaining mass to Mj-em = M — AM. 
Instead of requiring that the second member of the merger be a halo of mass M' = Mj-em, the 
Somerville & Kolatt algorithm allows for an arbitrary number of progenitors, each of which 
is chosen from Eq. ( 152|) . as well as an arbitrary amount of accreted mass. The next step is 
as follows. If Mrem > ^min, choosc another transition and get the new mass change AM'. If 
AM' > Minin, treat it as a merger and if AM' < M^i^, treat it as accretion. Update M^^ra 
by subtracting AM'. This procedure is then repeated until M^^m < ^min at which point 
the remaining mass is regarded as diffuse accretion. This produces a list of progenitors at 
a single timestep, though the list typically consists of two members and a small quantity 
of accreted mass. To fill out the tree, one applies this procedure with new timesteps to all 
the progenitors in the list and so on, recursively building a merger tree. The tree is of finite 
extent because each branch off of the tree terminates when all mass enters in units that are 
< Mmin, so that each unit is regarded as accreted material. 

The right panel of Figure [10] is reproduced from Figure 6 of Ref . . This panel shows 
the conditional mass functions of progenitors from an ensemble of merger trees generated 
using their algorithm compared to the analytic excursion set conditional mass functions. The 
suggestions of Somerville & Kolatt [2] have succeeded in drastically reducing the discrepancy 
between the analytic progenitor mass functions and the results of the Monte Carlo merger 
trees. 

Examples of mass accretion histories for main progenitors (the most massive progenitor at 
each timestep) computed according to this merger tree prescription are shown in Fig. [HI The 
large vertical jumps correspond to very sudden increases in mass due to (perhaps multiple) 
major mergers. The variety of different paths to the same final mass is evident. Also clear 
in Fig. [11] is the fact that smaller halos acquire their masses relatively earlier than larger 
halos, a picture that is consistent with the results shown in Fig. [7] Formation times are also 
easy to compute directly from these merger histories. The distribution of formation times 
among Monte Carlo realizations is shown in Fig. [12] Also shown in Fig. [12] are the formation 
time distributions from § IVII CI which appear to be in reasonable agreement (Note that the 



agreement here is better than that shown in the paper by Lacey & Cole [19| due to the use of 
the Somerville & Kolatt [2] algorithm for generating merger trees rather than a binary-split 
algorithm) . 

Merger trees also give a wealth of information about the manner in which the mass is 
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FIG. 11: Ten examples of mass accretion histories as a function of cosmological scale factor for 
main progenitors. These histories were generated using the Monte Carlo merger history method 
based on excursion set theory as described in the text. The solid lines correspond to a large cluster 
halo with a mass of 10^^ H^^Mq at z = (a = 1). The dashed lines correspond to a Milky 
Way-sized halo with mass M{a = 1) = 10-*^^ /i^^M©. 

acquired by the parent halos. Some properties of the distribution of all mergers are extremely 
simple. Consider the transition probability in the binary regime and in the regime AM <^ M 
so that the merging object is much smaller than the main progenitor with which it is merging. 
In this case, the transition probability is / ~ [Auj/{V27r\dS{M)/dM\)]AM-^/'^. As such, 
the average mass fraction that is accreted in units in units of mass in a logarithmic interval 
In(AM) varies as dF/dln(AM) ~ [Alu /{V27r\dS/dM\^/^)]AM^/^ for all timesteps. This 
requires only that the probability be in the binary regime and that the mass change be 
small. Integrating over all time, the full tree has this property as well with a deviation at 
AM ~ O.IM necessitated by mass conservation. The total mass accreted in units of mass 
AM per logarithmic interval dln(AM) integrated over the entire history of 1000 realizations 
of the merger history of an halo is shown in Fig.[T3l The scaling d-F/dln(AM) oc AM^/^ 
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FIG. 12: Halo formation times computed from the Monte Carlo merger tree algorithm compared 
to the closed form distribution of formation times given in § IVII CI The dashed lines correspond to 
halos of mass 10"*^^ H^^Mq, the solid lines correspond to halos of mass 10"*^^ H^^Mq = M^, and the 
dash-dot lines correspond to halos of mass 10^^ h~^MQ. The thin, smooth lines are the analytic 
predictions while the thick lines are the results 1000 realizations of the Monte Carlo procedure. 



at AM <^ M is evident. Integrating over In(AM) gives F = 1 because all mass is acquired 
in the form of halos of a particular mass in the excursion set model. 

Before moving on to extensions of the simple excursion set approach it is worth drawing 
attention to a simple but profound weakness of the excursion set theory that is not often 
discussed in the literature, though one noteworthy exception is the lecture by White jsoj]. In 
this and the preceding sections, I described how the excursion set theory is used to predict 
the merger rates between distinct dark matter halos. However, these distinct halos are never 
clearly identified in the excursion set approach. As an example, consider some patch of the 
initial Lagrangian space. We can draw a sphere about each point in that patch that has a 
smoothed density that is above threshold at some scale and thereby assign to each point a 
halo mass. In general, each point will be assigned a mass that is different from its neighboring 
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FIG. 13: The total fraction of a host halo mass at a = 1 accreted in units of size AM/M per 
logarithmic interval dln(AM), di<'/dln(AM). The result shown is from 1000 merger histories of a 
halo with mass 10^^ /j-^M© = il4 at a = 1. 



points. The result of this exercise is a continuous "halo mass field," rather than the definite 
assignment of mass in certain patches to distinct halos of particular masses. This makes 
tests of excursion set theory on an object-by-object basis a poorly-defined exercise because 
it is unclear how to assign individual particles to halos in the initial Lagrangian space. In 
fact, such tests have been made with different choices for assigning mass elements to halos 
(compare Ref. 59| and Ref. j4l|). This fundamental shortcoming of the excursion set theory 
may partially explain the difficulty of the theory in producing self-consistent merger trees 
and progenitor mass functions and its shortcomings compared to full numerical treatments. 



IX. BEYOND THE SIMPLE EXCURSION SET APPROACH 



In the proceeding sections, I have reviewed the simplest implementation of the excursion 
set approach to halo formation. While this approach provides valuable insight into the 
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influence of tlie statistical properties of the density field on the distribution and formation 
of halos, in its detail it is known to have a number of shortcomings, a number of which we 



have encountered already, that render it a 
is required (references {2, 2o|, 
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65, 



661] represent an incomplete 



sample of numerous studies that point out a variety issues). I will discuss briefly some of 
these shortcomings and some insights into the physical origin of these shortcomings in what 
follows. One point that is important to emphasize is that the excursion set approach is itself 
an approximation that neglects the details of nonlinear dynamics aside from the collapse 
threshold; however, the predictions discussed in the previous sections largely result from 
making several further assumptions that are motivated primarily by simplicity and may not 
be the basis of the most faithful implementation of excursion set theory. 



A. Nontrivial Barrier Shapes and Ellipsoidal Collapse 

One way in which the excursion set predictions are simplified is to assume that collapse 
and virialization of a halo occurs whenever the smoothed density exceeds a critical value 
6c, that is independent of any further details of the density field. Even some of the gross 
qualitative aspects of excursion set theory no longer hold when this simple form for the 
barrier is no longer required. 



Following Sheth 671], it is instructive to begin by studying a simple toy example such as the 
linear barrier where the threshold increases with the variance B{S) = 6c + f3S (constructed 
so that it can be easily solved). In this case the barrier recedes from the line 5 = at a 
"rate" (3. Solving for the first-crossing distribution of this barrier is equivalent to solving for 
the first-crossing distribution of a constant barrier where the transition probability causes 
trajectories to "drift" away from the barrier at a rate (3. In this case, we can revisit the 
logic that led to Eq. f l22l) but note that the mean transition should be (A5) = —f3AS so 
that Eq. fl23l) is replaced by 

dU _ 1 d'^U dU 

dS~2m^^'^^- ^^^^ 
The boundary conditions and the initial condition are the same as before. Note that in this 
case, the solution by using an image source of trajectories cannot be used because the image 
trajectories would not satisfy the same differential equation [Eq. ( |65l) ]. 

The solution to Eq. ( 165|) can be greatly simplified by making the transformation 11 = 
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U exp{P['j — PS/2]) where 7 = (5c — 5 as in §[IVl In that case U{S,'y) satisfies 

dU _ 1 d'^U 
'dS ~ 2~d^' 

which is the diffusion equation with no drift term, and the conditions 



(66) 



U{S = 0,7) = 5d(7 - 4) exp(-/57) (67) 

and 

?7(5,7 = 0) = (68) 

Following § HVt this gives the result 



exp — exp 



2S \ 2S 



(69) 



From Eq. ( I69|) . the fundamental relation for the probability for a mass element to be above 
threshold at a variance between S and S + dS is 

fUS\6o = 0,So = 0)dS = -i^ exp ( - ^-^^^) ■ (70) 



From this solution to the simple linear barrier a mass function can be derived in the usual 
way. 

It is useful to note a few things from this example. First, consider a mass function arising 
from Eq. (1701) . For S ^ Sc/(3, the mass function will have an exponential cut-off for high- 
mass (low-S) objects as before. Additionally, for 5* ^ 5c//5 the mass function will feature an 
exponential suppression of low-mass object as well [compared to f{S) ~ S^^^'^ for S >> 5^ 
in the constant barrier case, Eq. fl52]) ]. This is due to the fact that the barrier increases in 
height as S while the standard deviation of the walk increases only as so that the barrier 
becomes harder and harder to pierce at high S. In contrast, the constant barrier will always 
be crossed because the standard deviation of the walk v^, always becomes much larger than 
the barrier height as S increases. This leads to a second interesting feature. Integrating 
Eq. (1701) over all S gives the total fraction of mass in bound halos as F = exp(— 2/55c) < 1- 
The statement that all matter is bound into halos is not generally true in the excursion set 
approach, but is a result of specific barrier shapes. 

In fact, nontrivial barriers may play an important practical role in understanding things 
like the mass function and bias of halos. A number of papers have, essentially, argued 
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that the spherical collapse model that leads to the constant barrier at Sc is an inadequate 
description of the collapse of overdense patches and that the barrier shape can be more 
complex effectively to reflect both the more complicated collap se p rocess and the relative 
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differences between large and small overdense patches 
Tormen |34] noted that the GIF simulations obeyed a mass function where the fraction of 
mass in collapsed objects is modified from Eq. (!35|) to 

1 



Iv exp(— az/ /2), 



(71) 



with a = 0.707, p = 0.3, and A = 0.322 in order to guarantee that all mass is bound into 
halos of some mass. This is the form shown in Fig. [3] and it clearl y r epresents a major 



improvement over the simple excursion set mass function. Reference [3J] also showed that 
this form translated into a much improved model for halo bias (using the same logic as § |Vj) 
as shown in Fig. |H 



Subsequently, Sheth et al. [41] and Sheth & Tormen [42] discussed how the collapse 
fraction in Eq. (!7T|) is very similar to the collapse fraction one would derive from of an 
excursion set model with a nontrivial barrier, the form of which is motivated by considering 
the ellipsoidal, rather than spherical, collapse of overdense patches. In general, integrating 



the details of ellipsoidal col 



examples include Refs. 
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apse into the excursion set formalism is quite complicated (some 



72, iTSj). Sheth et al. [4l| greatly simplified this problem by 



considering the barrier that results form the mode of the distributions for ellipticity and 



prolateness of overdense patches (see Refs. 
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69| for details). In short, the simplification 



stems from assigning overdense patches only a single ellipsoidal shape rather than the full 
distribution of possible shapes that would be realized in the cosmological density contrast 
field. This yields a single barrier shape to describe ellipsoidal collapse rather than, at least 
effectively, considering a barrier with a shape that should depend upon the specific properties 
of each overdense patch. As it turns out, this approximation works fairly well. The resulting, 



simplified barrier is 
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cU 



5e(l+/30, (72) 

where the parameters f3 ~ 0.47 and 7 ~ 0.62 are determined by considering the evolution of 
an ellipsoidal density perturbation rather than a perfectly spherical density perturbation. 

Notice that the ellipsoidal collapse barrier of Eq. (!72|) is higher for low-mass (low-z/) 
objects. The physical reasoning for this is that low-mass objects typically exhibit greater 
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ellipticity {e.g., Refs. |4l|, |69j) and are more prone to disruption by tidal interactions so that 
their densities must be higher in order to overcome this countervaihng effect. Notice that 
one of the primary deficiencies of the constant-barrier excursion set mass function is that 
it predicts considerably too many low- mass halos (see Fig. [3]). As the ellipsoidal barrier is 
higher for higher 5* (lower mass), one should expect a reduction in low- mass halos relative 
to the standard excursion set mass function and this is precisely what is shown in Fig. [31 



Finally, Sheth & Tormen 42| showed that the ellipsoidal model still failed to model accu- 
rately conditional mass functions at small lookback times and large mass ratios. Moreover, 
these authors also pointed out that the model did a poor job with predictions of the number 
densities of halos in very overdense regions. This failure leads to the next simplification of 
the excursion set approach which I will address in § IIX CI However, before moving on, I will 
briefly note the methods used to solve the problem of a walk toward a generic barrier, as 
this problem seems to be increasingly interesting. 

B. General Barriers 

It can sometimes be useful to consider more general barrier shapes. For example, Furlan- 
etto et al. [t^ used a barrier with a more complicated shape to model the sizes of ionized 



regions during the epoch of reionization, while Sheth & Tormen [42| suggested that moving 
barriers could effectively encapsulate a wide variety of phenomena such as suppression of the 
collapse of small (low-mass) overdense patches in models where the dark matter is warm. In 
the latter case, the suppression is due to the non-negligible primordial velocity dispersions 
of the dark matter particles and the barrier shape could encode the importance of these 
velocities relative to the depths of the potentials. In this section, I describe two numerical 
methods to construct the first-crossing distribution for more general barrier shapes. 

The most commonly used method to solve for the first-crossing distribution of a barrier 
of more general shape is direct Monte Carlo construction of the distribution. This method 
is straightforward. One only needs to compute a large number of S{S) trajectories. To con- 
struct a trajectory, one computes values of 6{S) at a discrete set of values of 5*. Maintaining 
the assumption that the smoothing of the density fiuctuation field is performed with a k- 
space tophat filter, each step is independent and is drawn from Eq. ( l20l) for a step size AS". 
At this point, one only needs to keep track of the value of 5* at first crossing of the barrier. 
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FIG. 14: Monte Carlo solution for the first-crossing distribution of a linear barrier with /3 = 1, 
B{S) = 6c+S. The /e/t panel shows a sample of five S{S) trajectories that cross the linear threshold. 
The thin lines are the five trajectories and the thick line is the linear barrier. The histogram in the 
right panel shows the tabulated first-crossing distribution constructed from 10^ trajectories. The 
thin line shows the analytic solution to the first-crossing distribution from Eq. (j70p . It is clear that 
the Monte Carlo procedure reproduces the analytic result in the case of the linear barrier. 



As with generating halo merger trees, the choice of step size AS is important because all 
trajectories will be resolved only to within AS*. The step size must be small enough that 
a transition above the barrier is rare on any given step. This ensures that multiple barrier 
crossings are unlikely to be missed. 

I give an explicit example of the Monte Carlo method for the linear barrier with B{S) = 
6c + S in Fig. [TH The left panel shows a sample of five trajectories that cross the linear 
threshold. Of course, most trajectories do not cross the threshold in the range of S that I 
show. I have selected five that do and so these trajectories are a biased subsample of the 
ensemble of all trajectories. The right panel shows the tabulated first-crossing distribution 
constructed from 10^ trajectories. It is evident that the Monte Carlo method reproduces 
the analytic result of the previous section quite well. 

Recently, Zhang and Hui [3] described a very useful and elegant method for solving the 
first-crossing problem for barriers of a very wide range of possible shapes. Zhang & Hui 
pointed out that the problem can be recast in the form of an integral equation that can be 
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solved by well-known methods. 

The first step is to recognize that all trajectories either are included in 11(5, S) or penetrate 
the barrier so that 

/ f{S')dS'+ U{6,S)d6 = l, (73) 

Jo J-oo 

where f{S) is the first-crossing distribution and B{S) is the scale-dependent barrier height. 
If there were no barrier, the distribution of 11(5, S) would be given by \E'((5, S) from Eq. (!20l) . 
Therefore, 11(5, S) should be given by \E'(5, S) minus the trajectories that pierced the barrier 
at some S' < S and migrated downward to 5 by a variance S. Thus 

n(5, S) = ^(5, S) - f{S')^{S - B{S), S - S')dS'. (74) 
Jo 

Taking the derivative of Eq. (1731) with respect to 5* gives an integral relation for f{S) in 
terms of 11(5, S) and the barrier shape. Substituting 11(5, S) from Eq. ( 17il) into this integral 
relation for f{S) gives 

f(S) = G{S) + [' fiS')HiS, S')dS', (75) 
Jo 

where G{S) = [B{S)/S - 2dfi(S)/d5]^(5(5), 5), H{S,S') = [2dB{S)/dS - {B{S) - 
B{S'))/{S - S')]'^{B{S) - B{S'), S - S'), and H{S, S') = if 5' > S. 

Equation ( 1751) is an integral equation for the first-crossing distribution of f{S). Equa- 
tion ( |75l) has the variable as a limit of integration, and the unknown function of interest 
f{S) both under and outside the integral. This is the form of a Volterra Equation of the sec- 



ond kind. This equation can be solved numerically using standard methods {e.g., Ref. 76|). 
The first step is to discretize the integral on an interval from to S. The simplest approach 
is to choose a uniform step in S. Let n = 1,2, A^, AS* = S/N, Sn = nAS, fn = /(5'„), 
Gn = G{Sn), and Hnm = H{Sn, Sm — AS/2). Then, the integral can be rewritten as 

r^" . . AS ^ 

m=l 



r^" AS 

/ f{S')^{6 - B{Sn), S„ - S')dS' - — E Hnmlfm + fm-l], (76) 

JO 2 



so that 



. _ Gn + A^/2(Er=\ Hr,4fm + fm-i] + g(n,n)/„_i) 

^" l-iAS/2)H{n,n) ^ ' 



with the exceptions 



/o = Gq and 
r ^ Gi + (Ag/2)i7(2,l)[/i + /o] 
1- {AS/2)H{2,2) 
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FIG. 15: The first-crossing distribution for the quadratic barrier B{S) = 5c + 0.5S'^. The histogram 
is the result of a Monte Carlo construction of the distribution while the smooth, solid line is the 



iterative solution to the Volterra equation [Eq. ([75]) ] described by Zhang Sz Hui 



The values can be constructed iteratively based on Eq. (1771) . Utilizing Eq. (1771) as a 
solution for Eq. (175|) is generally a more efficient method to obtain a precise estimate for 
f{S) than computing a large ensemble of trajectories. 

As a simple example of this approach, I show the first-crossing distribution of the barrier 
B{S) = 6c + 0.55*^ in Figure HSl The histogram is the first-crossing distribution constructed 
from 10^ Monte Carlo trajectories, while the smooth line is the iterative solution to Eq. (175]) 
in Eq. (177]) . The two solutions are clearly in good agreement. Note also that the first- 
crossing distribution declines more rapidly at high S (small mass) than the first-crossing 
distribution for the linear barrier. This is a natural consequence of the fact that the barrier 
height grows more rapidly with respect to the variance than in the linear case so that the 
barrier is even more difficult to penetrate at high S. 
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FIG. 16: Correlated random walks. The thin, solid lines show five examples of correlated random 
walks of d{S) defined with a Gaussian filter in the left panel and a tophat filter in the right panel. 
These trajectories were constructed from five integrations of the Langevin equation [Eq. ()78p ] in 
each case. The thick line represents a threshold at 6 = 6c- The trajectories in this figure have been 
selected to pierce the barrier at S < 10 which is why they all tend to lie in the upper half plane. 
Notice that these trajectories vary smoothly, unlike the uncorrelated case. 



C. Correlations Between Scales 

As I described in § llVt the standard implementation of the excursion set theory makes 
use of a density contrast field filtered with a k-space tophat window function. This choice of 
a filter that is local in Fourier space requires that the mass that collapses onto a halo is not 
local in configuration space. In fact, the previous sentence understates the issue. In § HVl I 
noted that the use of the sharp k-space filter represents something of a problem because the 
sharp k-space window is difficult to associate with any particular mass because its volume 
diverges. However, once some method for assigning mass to a filter scale is decided upon, 
the use of the sharp k-space filter also has an important implication for the properties of 
halos. In the case of a sharp k-space window function, the transitions between different 
smoothed densities are independent Gaussian random variables as in Eq. ( |20ll . This implies 
that the formation histories of halos are completely unrelated to their local environments 
because in the excursion set formalism a halo of mass M has S{S) = 5c at S* = S{M) while 
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all steps away from this value are independent of each other. 

Though the effects are small, several studies have recently pointed out that the formation 
histories of halos in cosmological numerical simulations are correlated with their environ- 



ments 164 



791] as are internal properties of halos 
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79 



80l |. The implication of 



these studies is that correlations between scales cannot be neglected completely. Though 
the impact of these effects on observable quantities has not been thoroughly explored, the 



8l|. 



first observation of such an effect was recently presented by Berlind et al. 

Along other lines Sheth & Tormen 4^, attempted to determine whether or not their 
conditional mass functions were in error because of an incomplete treatment of ellipsoidal 
collapse and concluded that a more detailed treatment of ellipsoidal collapse could not ex- 
plain the failures of their model to predict conditional mass functions and density-dependent 
mass functions. They reasoned that because the worst failures occurred in high-density re- 
gions and for conditional mass functions at small lookback times and large mass ratios, the 
deficiency may be due to the neglect of correlations between different scales. The reasoning 
is as follows. Consider the specific example of the conditional mass function. To compute 
the conditional mass function, two smoothing scales are necessary. One scale describes the 
parent object while the other scale describes the progenitor halo at some earlier time. At 
high mass ratios and small lookback times, the two smoothing scales are quite similar so one 
might expect correlations to be important. However, in the low-mass-ratio regime where 
the progenitors are much smaller than the parent halo, the two smoothing scales are very 
different so that it might not be surprising for correlations between these very different scales 
to be relatively less important. 

Correlations between scales can be incorporated into the excursion set approach simply 
by modifying the window function. Any window function other than a sharp k-space window 
will result in walks with steps that are correlated to some degree. In the case of correlated 
walks, Eq. ( l2TI) no longer applies and so the results of § HV] are not valid. For correlated 
walks, the f{S)dS distribution must be computed by integrating numerically a Langevin 



18|. 



equation as described in the seminal work of Bond et al. in Ref. 

The Brownian motion of a particle can be treated by considering the acceleration of the 
particle as due to a fluctuating stochastic force F{t). The Langevin equation for the particle 
motion is dv/dt = F{t), where v is the particle velocity (for example, return to the review 
of Chandrasekhar in Ref. 23|]). In a similar way, the trajectory S{S) experiences a series 
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FIG. 17: The first-crossing distribution of the smoothed density field defined with a Gaussian 
window function. The left panel shows the first-crossing distribution explicitly. First the smooth, 
upper line is the analytic distribution for the sharp k-space window function given in Eq. p2p . 
The upper histogram shows the first-crossing distribution constructed by integrating the Langevin 
equation for this case and demonstrates that the two solutions are in good agreement. The lower 
histogram is the first-crossing distribution for the smoothed density field defined with a Gaussian 
window function. The right panel shows mass functions. The smooth line is the standard excur- 
sion set mass function. The lower histogram shows the mass function of halos using the Gaussian 
first-crossing distribution from the left panel and the S{M) relationship for a Gaussian filter func- 
tion. The upper histogram shows the mass function constructed from the Gaussian first-crossing 
distribution of the left panel but using the S{M) relation for a real-space tophat window function. 
The difference between the two histograms represents the effect of the different mass definitions 
with a fixed first-crossing distribution. 



of stochastic inputs from specific realizations of the different Fourier modes of the density 
field. The trajectory for the smoothed density field obeys 

^=g(lnfc)W^(fc;i?w), (78) 

where the variance S and the smoothing scale i?w are related through Eq. f|T^ as usual. 
The "stochastic force" in this case is Q{\n.k) which is a Gaussian random variable with zero 
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mean and a variance of 

(Q'(lnA:)) = ^ = A2(A:). (79) 

A large number of trajectories can be built up numerically by summing a large number 
of realizations of the stochastic forces. The algorithm is to select Q{lnk) on a grid of 
wavenumbers, to apply the smoothing to each Q{lnk) on a grid of smoothing scales, and to 
sum the contributions of each step in wavenumber to construct each trajectory. As before, 
the steps in smoothing scale should be sufficiently small that crossing the barrier of interest 
on any particular interval is unlikely, so it is often most efficient to use a non-uniform grid 
in resolution scale {e.g., see Ref. |18|). 

Figure [16] gives examples of trajectories with correlated steps computed by integrating 
numerically the Langevin equation of Eq. ( ITHl) . The left panel shows trajectories generated 
using a Gaussian window function and the right panel shows trajectories generated using a 
tophat window function. The most obvious feature of these trajectories is that they do not 
exhibit the frenetic variations of the uncorrelated random walks shown in Fig. [2] and Fig. [TH 
On the contrary, because all of the steps in the trajectory are correlated to some degree, the 
trajectories are smooth functions of scale. 

The first-crossing distribution for a set of trajectories with correlated steps can be con- 
structed numerically by generating a large number of trajectories like those in Figure [T6l 
Figure [T71 shows the first-crossing distribution for trajectories defined using a Gaussian win- 
dow function compared to the standard first-crossing distribution for sharp k-space filtering. 
In the range of S corresponding to halo masses of interest, the first crossing distribution of 
trajectories defined relative to a Gaussian window is lower than the standard, sharp k-space 
result. At large values of the variance {S^ 25), the first-crossing distribution of the trajec- 



tories defined relative to a Gaussian K 
features have been discussed in Ref. 



ter overtakes the sharp k-space distribution (these 
18|). Recall that translating f{S) into a mass function 
requires knowledge of the relationship between the variance S and the mass scale M. There- 
fore, comparing f{S) for two different filtering functions does not refiect the differences in 
the mass function directly because the two different S{M) relations must be accounted for. 
In the case of the sharp k-space filter, the definition of mass is far from obvious because 
the filter has a divergent volume associated with it. Standard practice is to use the S{M) 
relation set by a real-space tophat window function. The right panel of Figure [17] shows the 
mass function for the Gaussian first-crossing distribution using both the Gaussian window 
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relation for S{M) and the real-space tophat relation for S{M) so that the effect of changing 
the S{M) relation at fixed f{S) on the mass function is evident. 

D. Halo Formation Times and Environment: An Example 

The excursion set theory using filters other than the sharp k-space filter has been con- 
sidered only relatively rarely in the literature, but in light of recent results it is interesting 
to ask what excursion set theory might say about the relationship between halo properties, 
such as formation time, and halo environments. One can ask what predictions implementing 
a window that is localized in configuration space, rather than Fourier space, makes for the 
relation between, in particular, halo formation histories and large-scale halo environments. 

Consider as a concrete example the dependence of halo formation time on environment at 
fixed halo mass. A simple and physically reasonable way to introduce correlations between 
the past formation histories of halos and the environments in which halos reside is to consider 
trajectories of S{S) defined using a window function that is localized in configuration space. 
The logic for determining formation times is precisely the same as that in § IVII C\ but the 
probability distributions must be constructed numerically. I show the result of a simple 
toy example in Figure [181 The argument is fairly general, but for concreteness, take the 
specific example of the formation times of 10^^ h~^MQ halos with environment defined as the 
smoothed density field on scale of -Renv = 10 h~^Mpc, 6cnv = 5(-Rw = -Renv)- The calculation 
must be constructed numerically from a large ensemble of trajectories so it is necessary to 
take a bin of halo mass that has finite extent. For this exercise, I consider halos that make a 
first crossing of the barrier 6c between 1.3 < Si < 1.4. This corresponds to a mass range of 
approximately 1.8^ M/[10^'^ H^'^Mq]^ 2.2. This toy model is based on 10^ trajectories of 
S{S) defined using Gaussian filtering so that the steps of the random walks are correlated. 

The left panel of Figure [18] shows an explicit example of a conditional probability for 
crossing a second barrier of height u = 6c + Acj with Au; = 0.5 in an interval of width dS 
about S given a first crossing of the barrier at 6c in the range 1.3 < Si < 1.4, as specified 
in the previous paragraph. Also shown in the left panel of Figure [18] is the conditional 
probability for the quartiles of trajectories with the highest and lowest values of large-scale 
density 6env Notice that the conditional probability f{AS, Au!\Si,u!i) is shifted toward 
lower AS in regions of low density. Trajectories in low-density regions rise more rapidly 
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in the vicinity of the threshold. According to the logic of § IVII C\ the consequence of this 
is that trajectories in low-density regions generally will achieve higher values by the time 
the smoothing scale decreases to include half of the initial mass. As a direct result, these 
trajectories will generally correspond to halos with earlier formation times. I show this 
explicitly in the right panel of Figure [181 This panel shows the distribution of formation 
times computed from the same set of trajectories for all halos, as well as for the high and 
low quartiles of 5env 

Figure [18] shows clearly that introducing correlations between different scales into the 
excursion set formalism generally leads to halos in dense environments forming later and 
halos in diffuse environments forming earlier at fixed halo mass. The example of Figure [18] 
described in the previous paragraphs is fairly specific; however, while the magnitude of this 
effect varies, the sense of this effect is general and holds for different filter shapes and for 
(5env defined on widely different scales. Interestingly, I have computed the environment- 
dependent formation time distribution for halos of numerous masses, and the sense of this 
shift is independent of the mass of halo under consideration. 

The generality of this result implies that the excursion set theory with a density field 
smoothed by a filter that is localized in configuration space will lead to dense regions that 
preferentially host late-forming halos and vice versa for regions of low density. The result 
stems from the general fact that the conditional probability f{AS, Auj\Si,u!i) is shifted to- 
ward low AS* in low density regions. There is a simple qualitative explanation for this 
behavior. All trajectories are required to pass through the point Ui at Si in order to be 
considered a halo of a particular mass at all. Consider beginning a walk at S* = and 6 = 0. 
Trajectories that start so that the large-scale smoothed density is relatively low at some 
small value of the variance 5* must rise rapidly in order to cross Ui at 5*;. The steps in the 
walk or correlated, so that these trajectories typically continue to rise after crossing this 
threshold. The opposite is true in high-density regions. On average, trajectories that start 
at low large-scale densities, proceed to pierce subsequent thresholds earlier in the walk than 
counterpart trajectories with high large-scale densities. This general trend is illustrated in 
Figure [19] where 1 show both examples of trajectories in low- and high-density regions and 
a scatter plot of 5env against the rate of change of the trajectories about the crossing of 6c 
at Su d6{Si)/dS. 

The results of this toy model for the environmental dependence of halo formation indicate 
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that the natural expectation of excursion set theory should be that halos in dense environ- 
ments tend to form relatively later than average, while halos in diffuse environments tend 
to form relatively earlier than average. Excursion set theory with a constant barrier height 
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predicts this relative trend independent of halo mass. The studies in Refs. [77, 
indicate that low-mass halos (M -C M^) that form early should be found preferentially 
in high-density environments and conversely for late-forming halos. This is opposite the 
excursion set expectation described in the preceding paragraphs. On the contrary, the re- 
cent results of Wechsler et al. [78| and Wetzel et al. [79|] indicate that in the high-mass 



regime (M^ M^), it is indeed the late-forming halos that are preferentially located in dense 
environments. 

Wang et al. 8^ have recently studied the origin of the link between formation time and 
environment using a cosmological numerical simulation. These authors found that low-mass 
halos that form early are preferentially in dense environments, in agreement with previous 
work. More interestingly, Wang et al. continued to show that the accretion of matter onto 
small halos in dense environments is cut off prematurely (relative to excursion set expec- 
tations) due to tidal interactions and heating of the ambient environment by significantly 
more massive neighbors. The low-mass halos are forced to be early formers because their 
late accretion is eliminated in a competition with massive neighbors. Wang et al. propose 
this as the origin of the environment-dependent formation times measured in previous nu- 
merical studies. Finally, Ref. 82] went on to show that this competition-limited growth is 
not important for high-mass halos. 

It is interesting to conjecture that the prediction of the excursion set theory and the 
study by Wang et al. 8^ can be combined into a global picture that explains environment- 
dependent formation times at all mass scales. In such a scenario, the statistics of overdense 
patches in the primordial density contrast field would favor a bias toward late halo formation 
in dense regions. In the high-mass regime, halo growth would not be limited by competition 
in any meaningful way, so this bias would persist in the sense seen by Ref. [3] and Ref. [3] 
for high-mass halos in numerical simulations. However, in the low-mass regime, competition 
would stifle halo growth, forcing small halos in dense environments to be early formers. This 
gives rise to the bias of small halos in dense regions to be preferentially early formers. 

To my knowledge, the issues mentioned in this section have not been addressed in any 
adequate way at present. The model I have given here is an incomplete toy model intended 
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to illustrate the method for generalizing the excursion set predictions and to indicate the 
sense of the effect that excursion set theory predicts. A more detailed exploration of these 
points is an interesting undertaking, but is well beyond the scope of this review. 



E. The PINOCCHIO Algorithm 



In closing, I would like to point toward another tool for estimating halo mass functions and 

icly-available PINOCCHIO code 
16|'[85|]- (A similar tool that can be 
used to generate mock catalogs of halos rapidly is PTHalos developed by Scoccimarro and 



formation histories quickly and accurately. This is the p ub 



developed by P. Monaco and collaborators fly 



Sheth 
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8J|) The great value of the PINOCCHIO code is that it is considerably less compu- 
tationally expensive than a cosmological N-body simulation but predicts halo properties and 
formation histories that are in better agreement with N-body results relative to the simple 
excursion set approach (for a recent example see Ref. [661]). 

Briefly, the PINOCCHIO code uses an algorithm for identifying halos and halo properties 
such as mass, velocity, angular momentum, in particular realizations of the linear den- 
sity field. The PINOCCHIO code also contains a parameterized algorithm for building fila- 
mentary structures in the density field and for generating the merging histories of halos. 
PINOCCHIO creates a particle-based realization of the linear density field in a manner simi- 
lar to standard N-body simulations. PINOCCHIO then computes the approximate evolution 
of mass elements using a particular truncation of third order Lagrangian perturbation the- 
ory (for example, first order Lagrangian perturbation theory constitutes the well-known 
Zel'dovich Approximation used to initialize most N-body simulations). As PINOCCHIO is 
perturbative, the computation progresses at a fraction of the computational effort of an 
N-body simulation. PINOCCHIO groups mass elements together into halos using an orbit- 
crossing condition. When the orbits of two mass elements coincide, Lagrangian perturbation 
theory breaks down because the mapping from Eulerian position to Lagrangian position be- 
comes multi-valued and the density is infinite. At this point, the mass elements are grouped 
together into a halo (or perhaps a filament) in analogy with collapse to a point in the 
evolution of a spherical tophat overdensity. 

The development of the algorithm and the PINOCCHIO code involves a great deal of detail 
and is well-beyond the scope of this review. However, this development is well-documented 
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and the references above provide a good starting point for understanding the calculation. 
Acknowledgments 

I would like to thank Jeremy Tinker for making is numerical data available to me for 
this review. 1 thank James Bullock, Ncal Dalai, Anatoly Klypin, Audrey Kravtsov, Savvas 
Koushiappas, Brant Robertson, Eduardo Rozo, Douglas Rudd, Frank van den Bosch, Risa 
Wechsler, and David Weinberg for useful discussions and/or comments on the manuscript. 
This review is based on notes for lectures given at the Sixth Summer School of the Helmholtz 
Institute for Supercomputational Physics on Supercomputational Cosmology. The Summer 
School is funded by the Ministry for Science, Research, and Culture Brandenburg. 1 would 
like to thank the organizer of the Summer School on Supercomputational Cosmology, Ana- 
toly Klypin, for inviting me to give these lectures. I thank Anatoly Klypin, Audrey Kravtsov 
and Francisco Prada for encouraging me to expand them to the current review. I thank the 
Astrophysical Institute Potsdam for hospitality during my stay. This work was funded by 
the Kavli Institute for Cosmological Physics at the University of Chicago, by the National 
Science Foundation under grant NSF PHY 0114422, and by the National Science Founda- 
tion through the Astronomy and Astrophysics Postdoctoral Fellowship Program under grant 
AST 0602122. 



[1] W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974). 
[2] R. S. Somerville and T. S. Kolatt, MNRAS 305, 1 (1999). 
[3] S. Cole, C. G. Lacey, C. M. Baugh, and C. S. Prenk, MNRAS 319, 168 (2000). 
[4] R. S. Somerville and J. R. Primack, MNRAS 310, 1087 (1999), astro-ph/9802268. 
[5] G. KaufFmann and M. Haehnelt, MNRAS 311, 576 (2000), astro-ph/9906493. 
[6] E. Neistein, F. C. van den Bosch, and A. Dekel, ArXiv Astrophysics e-prints (2006), astro- 
ph/0605045. 

[7] A. R. Zentner and J. S. Bullock, ApJ 598, 49 (2003). 

[8] A. J. Benson, C. G. Lacey, C. S. Prenk, C. M. Baugh, and S. Cole, MNRAS 351, 1215 (2004), 
astro-ph/0307298. 

[9] S. M. Koushiappas, A. R. Zentner, and T. P. Walker, Phys. Rev. D 69, 043501 (2004), 



58 



astro-ph/0309464. 

[10] A. R. Zentner, A. Berlind, J. S. Bullock, A. Kravtsov, and R. H. Wechsler, Astrophys. J. , 
624, 505 (2005). 

[11] J. Penarmbia and A. J. Benson, MNRAS 364, 977 (2005), astro-ph/0412370. 
[12] J. E. Taylor and A. Babul, MNRAS 364, 515 (2005). 
[13] J. E. Taylor and A. Babul, MNRAS 364, 535 (2005), astro-pli/0410049. 
[14] R Monaco, T. Theuns, and G. Taffoni, MNRAS 331, 587 (2002). 

[15] P. Monaco, T. Theuns, G. TafFoni, F. Governato, T. Quinn, and J. Stadel, MNRAS 564, 8 
(2002). 

[16] G. Taffoni, P. Monaco, and T. Theuns, MNRAS 333, 623 (2002). 

[17] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986). 
[18] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Astrophys. J. 379, 440 (1991). 

[19] C. Lacey and S. Cole, MNRAS 262, 627 (1993). 

[20] C. Laccy and S. Cole, MNRAS 271, 676 (1994), astro-ph/9402069. 

[21] R. K. Sheth and G. Lemson, MNRAS 305, 946 (1999). 

[22] N. Wax, Selected Papers on Noise and Stochastic Processes (New York: Dover Publication, 

1954, edited by Wax, Nelson, 1954). 
[23] S. Chandrasckhar, Rev. Mod. Phys. 15, 2 (1943). 

[24] S. O. Rice, Bell Systems Tech. J., Volume 23, p. 282-332 23, 282 (1944). 

[25] S. O. Rice, Bell Systems Tech. J., Volume 24, p. 46-156 24, 46 (1945). 

[26] D. J. Eisenstein and W. Hu, Astrophys. J. 511, 5 (1999), astro-ph/9710252. 

[27] B. G. Williams, A. F. Heavens, J. A. Peacock, and S. F. Shandarin, MNRAS 250, 458 (1991). 

[28] S. M. Carroll, W. H. Press, and E. L. Turner, ARAA 30, 499 (1992). 

[29] S. Bildhauer, T. Buchert, and M. Kasai, A & A 263, 23 (1992). 

[30] R. A. Epstein, MNRAS 205, 207 (1983). 

[31] J. A. Peacock and A. F. Heavens, MNRAS 243, 133 (1990). 

[32] R. J. Bower, MNRAS 248, 332 (1991). 

[33] R. J. Adler, The Geometry of Random Fields (The Geometry of Random Fields, Chichester: 

Wiley, 1981, 1981). 
[34] R. K. Sheth and G. Tormen, MNRAS 308, 119 (1999). 

[35] A. Jenkins, C. S. Frenk, S. D. M. White, J. M. Colberg, S. Cole, A. E. Evrard, H. M. P. 



59 



Couchman, and N. Yoshida, MNRAS 321, 372 (2001). 

[36] N. Kaiser, ApJL 284, L9 (1984). 

[37] G. Efstathiou, C. S. Prenk, S. D. M. White, and M. Davis, MNRAS 235, 715 (1988). 

[38] H. J. Mo and S. D. M. White, MNRAS 282, 347 (1996). 

[39] U. Seljak and M. S. Warren, MNRAS 355, 129 (2004), astro-ph/0403698. 

[40] S. Cole and N. Kaiser, MNRAS 237, 1127 (1989). 

[41] R. K. Sheth, H. J. Mo, and G. Tormen, MNRAS 323, 1 (2001). 

[42] R. K. Sheth and G. Tormen, MNRAS 329, 61 (2002). 

[43] R. J. Schcrrcr and E. Bertschinger, Astrophys. J. 381, 349 (1991). 

[44] U. Seljak, MNRAS 318, 203 (2000). 

[45] C.-R Ma and J. N. Pry, Astrophys. J. 543, 503 (2000). 

[46] J. A. Peacock and R. E. Smith, MNRAS 318, 1144 (2000). 

[47] A. A. Berhnd and D. H. Weinberg, Astrophys. J. 575, 587 (2002). 

[48] V. Icke, MNRAS 206, 1 (1984). 

[49] E. Bertschinger, ApJS 58, 1 (1985). 

[50] R. van de Weygaert and E. Van Kampen, MNRAS 263, 481 (1993). 

[51] R. K. Sheth and R. van de Weygaert, MNRAS 350, 517 (2004). 

[52] G. R. Blumcnthal, L. Da Costa, D. S. Goldwirth, M. Lecar, and T. Piran, Astrophys. J. 388, 
234 (1992). 

[53] S. Furlanetto and T. Piran, MNRAS, in press (2005), arXiv:astro-ph/0509148. 

[54] S. Cole and N. Kaiser, MNRAS 233, 637 (1988). 

[55] S. Cole, Astrophys. J. 367, 45 (1991). 

[56] G. Kauffmann, S. D. M. White, and B. Guiderdoni, MNRAS 264, 201 (1993). 

[57] R. K. Sheth, MNRAS 281, 1277 (1996). 

[58] R. K. Sheth and J. Pitman, MNRAS 289, 66 (1997). 

[59] S. D. M. White, in Cosmology and Large Scale Structure, edited by R. SchaefFer, J. Silk, 

M. Spiro, and J. Zinn-Justin (1996), pp. 349-+. 

[60] J. M. Gclb and E. Bertschinger, Astrophys. J. 436, 467 (1994), astro-ph/9408028. 

[61] Y. P. Jing, ApJL 503, L9+ (1998), astro-ph/9805202. 

[62] G. Tormen, MNRAS 297, 648 (1998), astro-ph/9802290. 

[63] R. H. Wechsler, J. S. Bullock, J. R. Primack, A. V. Kravtsov, and A. Dekel, Astrophys. J. 



60 



568, 52 (2002). 

[64] R. K. Sheth and G. Tormen, MNRAS 350, 1385 (2004). 

[65] A. J. Benson, M. Kamionkowski, and S. H. Hassani, MNRAS 357, 847 (2005), astro- 
ph/0407136. 

[66] Y. Li, H. J. Mo, and F. C. van den Bosch, MNRASSubmitted, astro-ph/0510372 (2005), 

astro-ph/0510372. 
[67] R. K. Sheth, MNRAS 300, 1057 (1998). 
[68] P. Monaco, Astrophys. J. 447, 23 (1995), astro-ph/9406029. 
[69] J. R. Bond and S. T. Myers, ApJS 103, 1 (1996). 
[70] R Monaco, MNRAS 287, 753 (1997), astro-ph/9606027. 
[71] P. Monaco, MNRAS 290, 439 (1997), astro-ph/9606029. 
[72] E. Audit, R. Teyssier, and J.-M. Alimi, a 325, 540 (1997). 
[73] J. Lee and S. Shandarin, Astrophys. J. 500, 14 (1998). 

[74] S. Furlanetto, M. Zaldarriaga, and L. Hernquist, Astrophys. J. 613, 1 (2004). 
[75] J. Zhang and L. Hui, Astrophys. J. 641, 641 (2006), astro-ph/0508384. 
[76] G. B. Arfken and H. J. Weber, Materials and Manufacturing Processes (1995). 
[77] L. Gao, V. Springel, and S. D. M. White, MNRAS 363, L66 (2005). 

[78] R. H. Wechsler, A. R. Zentner, J. S. Bullock, A. V. Kravtsov, and B. Allgood, ApJ, In Press 

(astro-ph/0512416) (2006). 
[79] A. R. Wetzel, J. D. Cohn, M. White, D. E. Holz, and M. S. Warren, ArXiv Astrophysics 

e-prints (2006), astro-ph/0606699. 
[80] A. V. Maccio', A. A. Dutton, F. C. van den Bosch, B. Moore, D. Potter, and J. Stadel, 

MNRASSubmitted, astro-ph/0608157 (2006), astro-ph/0608157. 
[81] A. A. Berlind, E. Kazin, M. R. Blanton, S. Pueblas, R. Scoccimarro, and D. W. Hogg, ArXiv 

Astrophysics e-prints (2006), astro-ph/0610524. 
[82] H. Y. Wang, H. J. Mo, and Y. P. Jing, MNRAS, Submitted (astro-ph/0608690) (2006). 
[83] R. . Scoccimarro and R. K. Sheth, MNRAS 329, 629 (2002). 
[84] R. Scoccimarro, MNRAS 299, 1097 (1998). 



[85] The PINOCCHIO code can be obtained from URL http:/ /adlibitum.oats. inaf.it/monaco/Homepage/Pinocchio/ 



61 




FIG. 18: The dependence of halo formation times on halo environment in the excursion set theory. 
The figure shows results derived from trajectories of 6{S) defined using a Gaussian window function. 
The left panel shows the conditional probability for making a first crossing of a barrier of height 
uj = 5c + Auj with Alo = 0.5, given that the first crossing of the barrier Ui = 5c occurred in the 
interval 1.3 < Si < 1.4. This corresponds to a mass range 1.8 < M/llO^"^ H^^Mq] < 2.2. The solid 
histogram shows the full conditional probability distribution derived from 10^ such trajectories. 
The dotted histogram shows the conditional probability derived from the quartile of trajectories 
that have the highest smoothed densities on a scale of i?env = 10 /i~^Mpc. For the purposes of this 
exercise, the smoothed density on a scale of -Renv serves to define the environment. The dashed 
histogram shows the conditional probability from the quartile of trajectories with the smallest 
smoothed densities on a scale of i?cnv The right panel shows the distributions of formation times 
(in terms of the variable wp) derived from these trajectories. Again, the solid histogram shows 
the mean distribution of formation times for all halos, the dotted histogram shows the distribution 
of formation times for the quartile of halos with the highest large-scale densities, and the dashed 
histogram shows the formation time distribution for the quartile of halos in the lowest density 
environments. Notice that at fixed halo mass, halos with low large-scale densities typically form 
earlier (at higher uj-p) than halos that reside in regions that are overdense on large scales. 
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FIG. 19: Trajectories exhibiting low and high large-scale densities. The left panel shows twenty 
randomly-selected trajectories with values of 6env in the lowest quartile {solid lines) and 5cnv in 
the highest quartile (dashed lines). The horizontal, dotted line marks 5 = 5c indicating collapse at 
z = 0. All trajectories represent halos of nearly identical mass so they cross 6c at nearly identical 
points. The arrow indicates the variance on the length scale Rem/ = 10 /i^^Mpc which is where 
the environment measure (^env defined. Notice that trajectories with low large-scale densities 
tend to rise more rapidly after crossing the threshold at 6c- The right panel shows a scatter plot 
of large-scale density contrast 5env and the rate of change of 6{S) at the first-crossing of 5c at Si, 
d5{Si)/dS. The inverse correlation between large-scale density and local rate of increase is evident. 



